Symmetries and modelling functions for diffusion processes
aa r X i v : . [ phy s i c s . d a t a - a n ] J un Symmetries and modelling functions for diffusionprocesses
A.G. Nikitin a , S.V. Spichak a , Yu. S. Vedula b and A.G. Naumovets ba Institute of Mathematics of National Academy of Sciences of Ukraine,3 Tereshchenkivs’ka Street, Kyiv-4, Ukraine, 01601e-mail: [email protected]; b Institute of Physics of National Academy of Sciences of Ukraine,46 Prospect Nauki, Kyiv-28, Ukraine, 03028e-mail: [email protected]
Short title : Symmetries and modelling functions
PACS numbers : 68.43.Jk - Diffusion of adsorbates; 66.30.Dn - Theory of diffu-sion and ionic conduction in solids.
Abstract
A constructive approach to theory of diffusion processes is proposed, whichis based on application of both the symmetry analysis and method of mod-elling functions. An algorithm for construction of the modelling functions issuggested. This algorithm is based on the error functions expansion (ERFEX)of experimental concentration profiles. The high-accuracy analytical descrip-tion of the profiles provided by ERFEX approximation allows a convenientextraction of the concentration dependence of diffusivity from experimentaldata and prediction of the diffusion process. Our analysis is exemplified by itsemployment to experimental results obtained for surface diffusion of lithiumon the molybdenum (112) surface pre-covered with dysprosium. The ERFEXapproximation can be directly extended to many other diffusion systems.
Introduction
Experimental and theoretical studies of diffusion processes are of a great importancefor various branches of physics, biology, chemistry and other natural sciences. In ad-dition, such studies have important applications in medicine and many technologicalprocesses. A special interest is exited by surface diffusion processes which appear inmany physical and chemical systems. In particular, they are used in various kinds ofnanotechnologies.The theory of diffusion processes started in 1855 when Fick derived his classicaldiffusion equation [1] ∂θ∂t − ∂∂x a (cid:18) D ∂θ∂x a (cid:19) = 0 , (1)which still is a corner stone of the diffusion theory. In equation (1) D is a diffusioncoefficient, in general case depending on species concentration θ , and x a with a =1 , , a is imposed). Beingsupplemented by an appropriate initial data, equation (1) serves as a background fordescription of such diffusion processes which are characterized by diffusion flows linearin concentration gradients and not depending explicitly on space and time variables.Two standard problems of a diffusion theory are:1) To describe time evolution of the diffusion process, and2) To specify the dependence of the diffusion coefficient on concentrations of dif-fusing species.Of course, these problems are closely related, since if we know how the diffusioncoefficient depends on concentration θ , then the time evolution of the correspondingdiffusion process can be found using the Fick equation (1) and the related initial data.On the other hand, if we know θ as a function of time variable t and spatial variables x a , then we can find D solving the inverse diffusion problem using again equation (1).Both mentioned problems are very complicated and in general need rather sophis-ticated techniques. Even if we know the diffusion coefficient as an explicit function ofconcentration, then generally speaking it is possible to find only an approximate (nu-merical) solution of the first problem if at all. The second problem has a much morecomplex character, but in the case of a sharp step-like initial θ profile it is possibleto use the Boltzmann-Matano (BM) approach [2] and reconstruct the concentrationdependence D ( θ ) of the diffusion coefficient. This approach enables one to make anumerical calculation of the diffusion coefficient, but its accuracy is not very high,especially for small and large concentrations θ .Experimental data and numerical solutions are very important for description of adiffusion process, but to formulate its theory it is desirable to create some analyticalexpressions for studied values. Unfortunately, there are only few known exactly solv-able realistic diffusion problems, the most famous of them is probably the Barenblatone [3]. Thus it is a common practice to use rather rough analytic presentations of D ( θ ) to make a qualitative analysis of diffusion process (see, e.g., [4]).1n the present paper we propose a new method for description of time evolution ofa diffusion process and calculation of the diffusion coefficient. The distinct feature ofour approach is that we find both functions θ = θ ( t, x ) and D = D ( θ ) in an explicitform, i.e., solve both problems 1) and 2) analytically. To achieve this goal we startwith experimental data for a particular diffusion system and make the error func-tions expansion (ERFEX) of concentration profiles. Analytic description of diffusionprocesses is very convenient for their qualitative analysis. Moreover, our descriptionappears to be rather good quantitatively also; its deviation from experimental datadoes not exceed the inaccuracy of measurements.We apply this approach to describe in detail the surface diffusion of Li depositedon the molybdenum (112) surface which had been previously covered with a submono-layer of dysprosium. One more process, the diffusion of Dy adsorbed on Mo(112), isused to examine the method generality. Moreover, we believe that it can have a muchwider application area. Let us start with experimental data representing surface diffusion of Li on the Mo(112)surface precovered with a 0.25 monolayer of dysprosium (below we designate it asDy-Mo(112) surface). The data were obtained in ultra-high vacuum using local mea-surements of the work function by a contact potential method (see [5]-[7] for details).A schematic sketch of the method, termed scanning contact potential microscopy,is shown in Fig. 1. At the beginning of the experiment, we uniformly covered theFigure 1: Probing the surface distribution of adatoms by scanning contact potentialmicroscopy. (a) The initial step-like adatom distribution. (b) The adatom distribu-tion after surface diffusion.clean surface of a (112) oriented Mo single crystal with dysprosium (its surface densityamounted to 0.25 of a monoatomic layer) and equilibrated it by annealing at T=1100K. The Dy adatoms served thereafter as a controllable admixture which could affectthe diffusion kinetics of lithium [6]. Then, using a semiplane mask (screen) placedbetween the Li evaporator and the prepared substrate, a half of the crystal surface2as covered at room temperature with Li while its another half remained clean of Li(Fig. 1,a). The surface was scanned with an electron beam formed by an electron gun.The movable beam was used to record the distribution of the contact potential (workfunction) over the sample surface by P.A. Anderson’s (retarding field) method. Inseparate experiments, the work function change was carefully calibrated with respectto the absolute surface concentration of adsorbed atoms (adatoms). The calibrationdata served to convert the work function values to adatom concentrations, n . Tocharacterize the relative concentration of adatoms with respect to substrate surfaceatoms, we shall use, as it is conventional in surface science, the term ”degree ofcoverage” (or, for short, ”coverage”). It is defined as θ = n/n M , where n M is theconcentration of the substrate surface atoms ( n M = 8 . · cm − for the Mo(112)surface). The coverage θ = 1 is usually termed the geometrical monolayer.Under the experimental conditions provided in [5, 6, 7], there was neither evap-oration of the adsorbate into vacuum nor its diffusion (drain) into volume of theMo substrate. Thus the experimental results relate to a case a ”pure” surface diffu-sion which could be described by equation (1). Notice that the length of the crystalsample along the diffusion direction was about 10 mm while the extension of the dif-fusion zone in the experiments did not exceed 1.2-1.5 mm. Thus the boundary effectsconnected with the finite size of the sample could be neglectedThe initial Li distribution was step-like shaped. Then upon heating the adsorbateprofile spreads out due to surface diffusion. Since the edge of the step was orientednormally to the atomic channels on the Mo(112) surface, the diffusion proceededquasi-one-dimensionally along the channels, i.e. along the [1 , ,
1] direction [5-7].The experiment consisted of a series of measurements in which we recorded thetime evolution of the coverage profiles due to diffusion at a constant annealing tem-perature. At the beginning of each experiment (t=0), a standard (step-like) initialcoverage profile of the adsorbate was created and recorded on the crystal kept at roomtemperature, at which the adatom mobility is negligible. This profile is labeled t=0.Then the crystal sample was annealed at a fixed temperature. From time to time, theannealing was interrupted and the crystal quickly cooled down to room temperatureto record a new coverage profile arising due to diffusion. After that the annealing wascontinued, and so on. In this way we obtained a series of profiles corresponding todifferent annealing times and a constant annealing temperature. Such experimentscould be repeated at different annealing temperatures to determine the temperaturedependence of the diffusion kinetics.Measurements were made at times t = 0 , , , T = 600 K , experimental error in θ was ∆ θ = 0 . t = 0) is step-like shaped, but technologically it was impossibleto form an ideal step. The profiles obtained as a result of surface diffusion have acommon intersection point at x = 0 . , θ = 0 .
192 and have rather smooth contours.Moreover, a convention is used to set θ = 0 in the region where its value is below the3igure 2: Coverage profiles of Li adsorbed by Dy-Mo(112) at T=600 K: initial, t = 0(1), and measured at t = 1200 s (2), t = 2100 s (3), t = 3600 s (4), t = 5400 s (5).The x - coordinate gives distance in mm.measurement accuracy.Our task is to give a phenomenological theory of the related diffusion process.Abstracting from complicated underlying physical effects which are discussed in pa-per [8], we will describe the time evolution of the diffusion process and derive thedependence of the diffusion coefficient on concentration.To achieve our goal, we will exploit symmetries which form the integral part ofdiffusion processes. A precise analysis of the experimental data makes it possibleto find a specific symmetry which characterizes coverage profiles. Namely, let usfix a profile θ recorded at time t and consider it as a given function of x . Thenany other profile θ ′ measured at time t ′ can be obtained from θ using the followingtransformation: x → x ′ = ( x − x ) r tt ′ + x , (2)where x = 0 .
88 is coordinate of the point common for all coverage profiles.This statement can be verified directly or using computer fits to the experimentalcurves. For example, starting with experimental data for t = 3600 s and applyingtransformation (2), one can reproduce profiles for t = 1200 s , s and 5400 s .It should be stressed that the found symmetry is rather exact. As a rule, adeviation of curves obtained using transformation (2) from experimentally measuredprofiles is within the limits of experimental error, and this deviation decreases withgrowing time. For example, comparing experimental data for coverage profile at4 = 5400 s (Fig. 2) and the corresponding values obtained via change (2) we concludethat they are very similar, namely, the differences are below the experimental error. The exact meaning of the symmetries discussed in Section 2 is that functions θ ( t, x )describing coverage profiles are invariant with respect to the following one-parametricgroup of transformations: t → t ′ = e α t, x → x ′ = xe α + x (1 − e α ) , (3)where α is a real parameter. Indeed, solving the first of equations (3) for e α and usingthe second equation we come to relations (2).Starting with (3) and using tools of the classical group analysis, it is possibleto describe time evolution of profile θ ( t, x ). Indeed, the infinitesimal operator oftransformations (3) has the form: X = η ∂∂t + ξ ∂∂x ≡ t ∂∂t + x ∂∂x − x ∂∂x , where η = ∂t ′ ∂α | α =0 and ξ = ∂x ′ ∂α | α =0 [9]. The invariance of coverage profiles with respectto transformations (3) means that θ ( t, x ) solves the following equation [9]: Xθ ( t, x ) = 0 . (4)It follows from (4) that time evolution of coverage profiles θ ( t, x ) is described bythe following equation: ∂θ∂t = x − x t ∂θ∂x or ∂ ˜ θ∂t = y t ∂ ˜ θ∂y , (5)where y = x − x , ˜ θ ( t, y ) = θ ( t, y + x ). As an initial condition we can choose one ofmeasured profiles, say that one which corresponds to t = t = 1200:˜ θ ( t , y ) = θ = θ ( y + x ) (6)where θ ( x ) is a function given numerically in Table 1. Henceforth we omit tilde andwrite θ ( t, y ) instead of ˜ θ ( t, y ).Thus we can describe the time evolution of the coverage profiles without a diffusionequation. Solving problem (5) with condition (6) and using numerical data given inthe Appendix we can find shapes of these profiles for any time. Some of such profilesare presented in Fig. 3.As we see, the very existence of the symmetries in experimental data makes itpossible to predict shapes of the coverage profiles which will appear at various times5igure 3: Coverage profiles θ ( t, x ) for Li on Dy-Mo(112) at T=600 K: experimentalfor t = 2100 s (1) and theoretical for t = 7000 s (2), 10000 s (3), 20000 s (4), 100000s (5). The x -coordinate gives distance in mm.of heating. This statement is valid for any diffusion system which admits symmetries(3). But if we are interested in analytical description of the diffusion process, wehave to pose initial conditions analytically. A basic problem is to create a consistentmodel of the diffusion process, i.e. find the dependence of the diffusion coefficient onconcentration, which is important for physical interpretation.In the following sections we solve this problem and give explicit representationsof coverage profiles θ ( t, x ) in analytical form. Transformations (2) have a stable point x = x = 0 .
88 and θ = 0 .
194 which was theonly one not being changed. It is convenient to choose just x as a zero point of ourcoordinate system, i.e. to use variable y = x − x = x − .
88 instead of x .Taking into account that in fact we deal with a process which is one-dimensionalwith respect to spatial variables, it is possible to reduce equation (1) to the followingform: ∂θ∂t − ∂∂y (cid:18) D ∂θ∂y (cid:19) = 0 (7)where D = D ( θ ) , y = x − x .Formula (7) represents a rather complicated non-linear partial differential equa-tion. In addition, we do not know the θ -dependence of the diffusion coefficient D .6ortunately, this equation has a very useful symmetry with respect to scaling of in-dependent variables, being invariant with respect to transformations (3), or t → t ′ = e α t, y → y ′ = ye α , θ → θ ′ = θ. (8)Just this nice property of the diffusion equation makes it possible to use theBoltzmann variable ξ = y √ t and search for its similarity solutions θ = θ ( ξ ) whereboth θ and ξ are invariants of transformations (8).Notice that equation (7) is invariant also with respect to shifts of independentvariables t → t + b , y → y + k (9)with arbitrary real parameters b and k . This invariance allows one to choose arbitraryinitial time and justifies the transition from x to y which we made above. For someparticular functions ˜ D ( θ ), symmetry of the Fick equation is more extended. A com-plete classification of symmetry groups of equation (7) has been made by Ovsiannikov[10], a complete group classification of systems of two diffusion equations with sourceterms can be found in [11] and [12].However, symmetries (8) should be compatible also with the initial data of ourproblem, which are of the form: θ (0 , y ) = θ ( y ) , (10)where θ is the initial coverage profile ( t = t = 0) represented numerically in theAppendix and graphically in Fig. 2.We see that the experimentally created profile at t = 0 is not strictly step-shaped,and so is not invariant with respect to scaling (8) . However, the profile θ is not farfrom a step and can be considered as a perturbed Heaviside function H ( − y ) multipliedby 0 .
327 ( θ max = 0 .
327 is the maximum coverage in the initial θ profile). On the otherhand, and it is an experimental fact, for sufficiently large times t solutions of ourproblem indeed are invariant with respect to transformations (2). And if we applythis transformation to infinitely small t ′ , all profiles θ , θ ...θ tend to a step-shapedone. So we have a direct experimental confirmation of the known mathematical factthat similarity solutions of the diffusion equation (7) can serve as attractors for othersolutions. In other words, the yet unknown diffusion coefficient should depend on θ insuch a way that the related Cauchy problem (7), (10) be asymptotically stable withrespect to small perturbations of the initial data (see [13] for exact definitions).Thus, instead of the actual initial data presented in the Appendix and Fig. 2, wecan consider an idealized situation when the initial coverage has a step shape, andto suppose that for t = 0 the concentration is proportional to the Heaviside function H ( − y ): θ (0 , y ) = θ max · H ( − y ) = (cid:26) . , y < , , y ≥ θ max = 0 .
327 is the maximum concentration in the initial coverage profile.The initial-value problem (7), (11) is invariant with respect to transformations (8),and so it is possible to search for invariant solutions θ = θ ( ξ ) depending on invariantvariable ξ = y √ t . In this way the problem is reduced to the following one: ξ ∂θ∂ξ + 2 ∂∂ξ (cid:18) D ∂θ∂ξ (cid:19) = 0 ,θ ( −∞ ) = 0 . , θ (+ ∞ ) = 0 . (12)Remember that we do not know yet the dependence of the diffusion coefficient D on degree of coverage θ . Using experimental data describing dependence of θ on x atfixed heating time t and applying the BM approach [2] it is possible to calculate D numerically. Unfortunately, such calculations cannot be done with a sufficiently goodaccuracy, especially in the region of small concentrations. In addition, in this way wecannot find the diffusion coefficient in an analytical form. In Section 6 we suggestanother approach which presupposes direct analytical modelling of profiles θ ( ξ ). Thus we had formulated a possible model of the analyzed diffusion process which isbased on equation (7) whose solutions should satisfy the idealized initial conditions(11). However, the possibility of using Fick equation (7) is nothing but a suppositionwhich needs additional justifications. In particular, it is necessary to be ensured thatthe diffusion flow is linear in the concentration gradient.However, we can use the well justified fact, i.e., the invariance of experimentalcoverage profiles with respect to transformations (2), and set the following problem:to find the most general evolution equation which is compatible with this symmetry.Thus let us suppose that the evolution equation admits symmetries (2) and alsoshifts (9) (i.e., does not depend explicitly on space and time variables). Then usingtools of classical Lie analysis [9], we easily find its general form: ∂P∂t − ∂∂y (cid:18) ˜ D ∂P∂y (cid:19) + G (cid:18) ∂P∂y (cid:19) = 0 (13)where P is a dependent variable whose evolution we need to describe, ˜ D and G arearbitrary functions of P . In particular P can represent a degree of coverage, in thiscase it should be changed to θ .We do not present the related routine calculations since our rather strong re-strictions (equation should be of evolutionary type and admit the above mentionedsymmetries) in fact reduce the problem of deducing of (13) to direct use of the di-mension analysis.Equation (13) is a direct generalization of the Fick equation (7), and the lattercorresponds to a particular choice G = 0. 8hus starting with symmetries (2), which can be found in the experimental data,and symmetries (9), which are natural transformations for the considered diffusionsystem, we deduce the generalized diffusion equation (13) which includes the Fickequation (7) as a particular case. Possible physical motivations for generalization ofequation (7) to (13) are discussed in Section 10. It is well known that in the case when diffusion coefficient D is independent onconcentration, the general solution of the problem (7), (11) is given by the followingequation: θ ( t, x ) = a erfc ( bξ )where a = θ max / , b = 1 / (2 √ D ), ξ = x/ √ t and erfc is the complementary errorfunction: erfc ( z ) = 1 − erf ( z ) , erf ( z ) = 2 √ π Z z e − t dt. This fact suggests using just the error function as a constructive element of concen-tration profiles for θ -dependent diffusivities. This idea appears to be very successfulfor description of diffusion processes in general and the processes discussed in Sections2 and 9 in particular.In this section we present and discuss examples of modelling functions for the cov-erage profiles. An algorithm for constructing such functions is given in the followingsection.First we consider a rather straightforward representation of profiles θ ( ξ ) which,however, is valid only in a reduced interval of the Boltzmann variable ξ : θ = θ (1) = a (1 − erf ( b ξ ) ) (14)where a and b are parameters. Asking for minimal mean-square deviation of function(14) and using MAPLE tools we fix parameters a and b to be a = 0 . , b = 0 . . Function θ (1) perfectly describes the shape of profile θ ( ξ ) for ξ lying in the interval0 . < ξ < ∞ , see Fig. 4. In this interval the discrepancy | θ − θ (1) | does not exceedinaccuracy of measurements. However, for ξ < . ξ , it is sufficient toadd a small extra term to θ (1) and define: θ (2) = θ (1) + a erfc ( b ξ ) , ξ ≥ θ ( ξ ) for t = 5400 (full curve) and the curve θ givenby relation (14) (broken curve). Units for ξ are 10 − (mm/s / )where a = 0 .
02 and b = 1 . t = 3600 s; the deviation | θ − θ (1) | does not exceed the inaccuracy ofmeasurements. Fig. 5 presents the experimental curve and curve defined by equation(15), and it is seen that they practically coincide. More exactly, the deviation ofvalues of function (15) from experimental data is less than 0 . θ ( ξ ): experimental data for t = 3600 s (full curve) andfunction given by relation (15) (broken curve). Units for ξ are 10 − (mm/s / )One more possible modelling function is given by the following equation: θ (3) = a ′ (1 − erf ( b ′ ξ ) ) + a ′ erfc ( b ′ ξ ) (16)where a ′ = 0 . , b ′ = 0 . , a ′ = 0 . , b ′ = 0 . t = 5400 s is less than 0 . a priori supposition that the coverage profiles can be described by second orthird order polynomials of error functions and demanding minimal root-mean-squaredeviation of these polynomials values from the coverage profiles. This supposition isnot necessarily valid for other diffusion systems, e.g., it does not give a constructiveway to build approximate modelling profiles for the system discussed in Section 9.In the following section we present a regular way to calculate modelling functionsfor any diffusion system using the error function expansion (ERFEX). We have shown above that erfc functions can be successfully used as construction ele-ments of the modelling functions of coverage profiles at least for a particular diffusionprocess. Now we shall give a regular way for constructing such functions which has amuch more extended application area.Of course a concrete form of the modelling functions strongly depends on thediffusion system, and we cannot propose a universal method how to obtain the mostsimple and exact analytical form of an arbitrary concentration curve. Nevertheless,in this section we give an algorithm for calculation of the modelling functions whichcan be applied to any sufficiently smooth concentration profile in the diffusion zone.In general case we propose to use ERFEX and search for modelling functions inthe form θ ( ξ ) = n X i =1 A i erfc ( k i ( ξ − ξ i )) (17)where ξ i ( i = 1 , , ...n ) are some fixed values of Boltzmann coordinate, k i and A i areparameters which should be specified.Let ξ i − < ξ i , or ξ < ξ < · · · < ξ n , then optimal values of k i lie inside the interval0 . ξ i +1 − ξ i ≤ k i ≤ ξ i +1 − ξ i . (18)In particular it is possible to choose the points ξ , ξ · · · , ξ n in a regular way, i.e.,with the same distances ξ i +1 − ξ i for all i , and to fix all parameters to be k = k = · · · = k n = p with some p compatible with (18).Let θ , θ ...θ n be known values of coverage at points ξ , ξ , ..., ξ n and M be a matrixwhose elements are M ij = erfc ( k i ( ξ i − ξ j )). Then parameters A , A , ..., A n are easilyfound by solving the following system of linear algebraic equations: M ij A j = θ i , i = 1 , , · · · n, (19)11here summing up over the repeated index j is imposed from j = 1 to j = n .Equations (19) are nothing but relation (17) considered at points ξ , ξ , ..., ξ n .Let us apply the algorithm to find a modelling function for the diffusion pro-cess discussed in Sections 2-6. Consider the results presented in the Appendix for t = 5400 s. For simplicity we choose points ξ , ξ , .. in a regular way. Namely, wechose in the table selected data corresponding to x = 0 .
63, 0 .
68, 0 .
73, 0 .
78, 0 . .
88, 0 .
93, 0 .
98, 1 .
03 and 1 .
118 (remember that x = 0 . k i = 0 . i = 1 , , ...n ( n = 12). At that, system (19) is easily solved with using symboliccalculus program MAPLE. As a result we find the following representation for θ ( ξ ): θ ( ξ ) = − . erfc (0 . ξ + 3 . . erfc (0 . ξ + 3 . − . erfc (0 . ξ + 2 . . erfc (0 . ξ + 1 . . erfc (0 . x + 1 . . erfc (0 . x + 0 . . erfc (0 . ξ − . . erfc (0 . ξ − . . erfc (0 . ξ − . . erfc (0 . x − . − . erfc (0 . ξ − . . erfc (0 . ξ − . . (20)This a bit cumbersome formula appears to be very precise; in the interval − <ξ < . t , it is possible to find ”the most exact” modelling functions θ E whose valuessimply coincide with the experimental curve. However, taking into account the valueof experimental error, such business does not seem to be reasonable.Notice that using the algorithm with a non-regular distribution of points ξ i it ispossible to find a more simple form of the modelling function: θ ( ξ ) = X i =1 A i erfc ( Rξ + B i ) , (21)where A = 0 . , A = 0 . , A = 0 . , A = 0 . ,B = − . , B = − . , B = 1 . , B = 3 , R = 0 . . (22)Being much more simple than (20), function (21) is rather exact too; the corre-sponding standard quadratic deviation from experimental results is less than 0 . Thus we have at our disposal analytic expressions for a close approximation of cov-erage profiles. Using them we can calculate the diffusion coefficient D ( θ ) by directintegration of equation (12). 12irst using the simplest modelling function (14) we find D ( θ ) for concentrations θ < .
16, the related values of ξ satisfy ξ > .
6. Substituting (14) into (12) andintegrating the resultant expression from ξ = z > a √ πb (cid:18) erf ( b ξ ) e − ( b ξ ) + 1 √ erfc ( √ b ξ ) (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) z ∞ − (cid:18) a b √ π erf ( b ξ ) e − ( b ξ ) D (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) z ∞ = 2 a √ πb (cid:16) erf ( b z ) e − ( b z ) + 4 b erf ( b z ) e − ( b z ) D (cid:17) = 0 . (23)Solving (23) for D and using (14) we obtain: D = 14 b e Z erfc ( √ Z ) q − θa ) , Z = erfinv r − θa ! . (24)where erfinv ( · ) is the inverse error function defined by means of θ = erf ( erfinv ( θ )) . Formula (24) defines the diffusion coefficient D as an explicit function of concen-tration and is valid for all θ lying in interval 0 ≤ θ < . D starting with other modelling functions θ ( ξ ) found inSections 6 and 7. The general expression for D is given by the following equation: D ( ξ ) = R ∞ ξ z dθ ( z ) dz dz dθ ( ξ ) dξ , (25)which, together with a modelling function for θ ( ξ ), determines the diffusion coefficientas a function of θ given in a parametric form.In contrast to the BM approach, to find the dependence of the diffusion coefficienton concentration we simple need to calculate a definite integral of the known functionand divide it by the redoubled derivative of (known) function θ with respect to ξ . Allthese operations are easily handled using, e.g., MAPLE tools.In Fig. 9 we present a plot of the diffusion coefficient (25) versus the degree ofcoverage θ given by relation (20). These D values are consistent to within 10-25%with the data obtained in work [6]. However, the maximum at θ ≈ .
25 was notrevealed in [6], where the graphical Matano’s evaluation of the coverage profiles wasapplied. This demonstrates that the ERFEX approach provides a more accurateprocessing of experimental diffusion results.13 ,05 0,10 0,15 0,20 0,25 0,302345678910
Figure 6: Diffusion coefficient D ( θ ) ( in units mm / s · − ) of Li on Dy-Mo(112)calculated using modelling function (20). T=600 K. We have seen that modelling functions (20), (21) give very exact analytic expressionsfor coverage profiles of Li adsorbed on Mo(112). A natural question arises whether itis possible to find such functions for modelling the profile shapes for other systems.In this section we apply the ERFEX method to another adsorption system, namelyto Dy adsorbed on Mo(112). Experimental data for this system were obtained anddiscussed in [7]. The related plots of coverage profiles are presented in Fig. 7.We see that all the profiles recorded in the diffusion process again have a commonintersection point this time at x = 1 .
69 which, however, lies out of the initial profile.Their shapes are much more specific than ones given in Fig. 2. There are two zoneswith a quick change of coverage and three zones where this change is rather slow.These profiles mirror a structural self-organization in the diffusion region, i.e. forma-tion of a series of two-dimensional adsorbate phases which differ from each other bydiffusion parameters and mechanisms [7]. Nevertheless, it appears possible to describethese profiles analytically.First we represent the coverage profiles using the Boltzmann variable ξ = x/ √ t and setting the reference frame x = 0. As a result we conclude that profiles for t = 2400 s and t = 4800 s became rather close, the mean quadratic deviation of themis less than 0.003. Thus it is possible to describe time evolution of profiles using theapproach discussed in Section 3. 14igure 7: Coverage profile of Dy adsorbed by Mo(112) at T=800 K: initial, t = 0 (1)and measured at t = 360 s (2), t = 2400 s (3), t = 4800 s (4) [7].Using again the ERFEX, we find the following representation for the coverageprofile measured at t = 4800 s: θ ( ξ ) = 0 . erfc (2 . ξ − . . erfc (0 . ξ − . . erfc (0 . ξ − . . erfc (4 . ξ − . . erfc (2 . ξ − . . erfc (1 . ξ ) + 0 . erfc (0 . ξ +2 . . erfc (0 . ξ + 5 . . erfc (1 . ξ + 15 . . erfc (2 . ξ + 22 .
2) + 0 . erfc (1 . ξ + 19 . . (26)Function (26) is a rather precise approximation of the coverage profile θ ( ξ ) at t=4800s; the mean quadratic deviation from the experimental results is less than 0.0015. Aplot of the related curves is given by Fig. 8.We see that ERFEX provides a very good analytical representation for a compli-cated coverage profile of Dy adsorbed on Mo(112).
10 Discussion
The theory of diffusion is both very old and good developed [14]. Nevertheless, it stillcontains a lot of unsolved problems which attract attention of numerous investigators.In the present paper we study three aspects of this theory: using of symmetriesin experimental data to describe the time evolution of a diffusion process, search-ing for a generalized Fick equation which is compatible with these symmetries, andconstruction of modelling functions to describe concentrations of diffusing substancesand calculate the diffusion coefficient.Like the diffusion theory, the basic branch of mathematics which deals with sym-metries, i.e., the theory of continuous groups, is rather old too. It was started by15igure 8: Coverage profile of Dy adsorbed on Mo(112) versus the Boltzmann variable,measured at t = 4800 s (full curve) and profile described by function (26) (dottedcurve). T=800 K.Sophus Lie around 20 years after appearance of the Fick theory. The classical group analysis and its modern generalizations present effective toolsfor investigation and applications of symmetries of mathematical models, includingdiffusion ones. However, to apply these tools directly it is necessary to start with amodel formulated in terms of (partial or ordinary) differential equations. A specificity of diffusion systems is that in general we do not know evolution equa-tion a priory. Even if there is a cogent argumentation that the process is described bythe Fick equation, the dependence of diffusion coefficient on concentration is usuallyunknown. Thus to apply the Lie analysis it is reasonable to search for symmetries inexperimental data . And it is the first idea which we use in the present paper.The second idea is to use erfc functions as a constructive elements of modellingfunctions for concentration curves, i.e., to apply the ERFEX expansion. There aretwo origins of this idea: first, a similarity to solutions of the Fick equation witha constant diffusion coefficient, where the erfc function appears very naturally, andsecondly, the specific shape of this function which seems to be ideal for its use as a”brick” for building typical concentration curves.Our research is based on a particular diffusion system, i.e., Li adsorbed on Dy-Mo(112), which is well studied and described in papers [5]-[8]. This system represents The background of the continuous group theory was formulated by S. Lie in 1875 and publishedin 1876 [15] Integral, difference and fractional differential equations also are subjects of modern group anal-ysis θ andthe relevant system property P are related by the dependence θ = F ( P ) . (27)Substituting (27) into the Fick equation (7) we come to equation (13), where˜ D ( P ) = D ( F ( P )) , G = ˜ D F ′′ F ′ , F ′ = ∂F∂P . (28)On the other hand, if a diffusion process is accompanied by dissipation (e.g. due toevaporation or chemical reaction of the diffusing species), the related mathematicalmodel also needs a generalization of the Fick equation by inclusion of the termsdepending on ∂θ∂x . In order to create a mathematical model of such process whichkeeps symmetries (3), one should use just the motion equation (13) with P = θ .We see that the generalized equation (13) appears as a result of rather straightfor-ward considerations and presents an alternative to the generally recognized equation(7). Surely, for the case when θ is a linear function of P and there is no dissipation,equations (13), (28) are reduced to (7).Use of modelling function for coverage profiles is very convenient for analysis ofdiffusion processes. These functions are sufficiently exact and can be used for directcalculation of diffusion coefficients starting with the Fick equation.In the present paper we propose a few modelling functions for coverage profilesof Li adsorbed on Dy-Mo(112). They described experimental data with a high preci-sion, the deviation from experimental data is less than experimental error. Thus themodelling functions present useful tools for both qualitative and quantitative studiesof the diffusion systems.Let us note that calculating the modelling function for small concentrations withhigher accuracy (i.e., using all experimental given points) it is possible to find a spikefor the diffusion coefficient for small concentrations. The related plot is given by Fig.9. 17 ,000 0,001 0,002 0,003 0,004 0,005 0,006 0,007 0,0080,51,01,52,02,53,03,5 Figure 9: Diffusion coefficient D ( θ ) ( in units mm / s · − ) of Li on Dy-Mo(112)calculated precisely for small concentrations. T=600 K.Notice that the spike in Fig. 9 is indicated for a very small concentrations com-patible with the experimental error and so it cannot be treated as well experimentallyjustified. On the other hand this spike can be observed for all series of experimentaldata presented by Fig. 2 and also in some other diffusion systems, e.g., in Dy ab-sorbed on Mo(112). Indeed the coverage profile given by Fig. 8 is rather steady at ξ = 5 − θ ∼ .
02 and so the related diffusion coefficient will have a maximumthanks to small value of the derivative ∂θ/∂ξ .Concerning interpretation the spike indicated by Figs. 9 and min.-max. presentedby Fig. 6 we can mention that the coverage dependence of diffusivity is due to lateralinteractions of diffusing atoms and partially also to the specific features of the sub-strate atomic structure (both intrinsic and caused by defects). The combined actionof lateral interactions and surface potential corrugation determines the sequence ofphase transitions that occur in the diffusion zone. The phases can differ from one an-other not only by the diffusion parameters, but also by the diffusion mechanisms (seee.g. a review [8] and references therein). The maxima of diffusivity at low coveragewere observed for a number of systems. We attribute this effect to a high mobilityof single adatoms. The fast decrease of D with growing coverage may be causedby formation of clusters, whose diffusion mechanisms can be very diversified, buttheir mobility is generally lower than that of single adatoms. The diffusivity is alsorather low in the regions of first-order phase transitions in adlayers. As the coveragegrows further and the adlayer becomes increasingly dense approaching a close-packedmonolayer, the diffusion takes on a pronounced collective character. In particular,18n the region of commensurate-incommensurate phase transition the diffusion seemsto proceed by a relay-race walks of misfit dislocations (topological solitons), whichprovides a high diffusion rate (a local D maximum). This effect was revealed for manyadlayers. For more details refer to papers [8], [16].Summarizing, we propose a constructive and convenient algorithm (ERFEX) forgenerating of modelling functions which is valid for arbitrary sufficiently smoothcurves not necessarily related to a diffusion process. In particular, using this al-gorithm and starting with experimental data, it is possible to determine the diffusioncoefficient with a higher accuracy than the BM approach and spline approximation.The algorithm can be treated as a specific generalization of the wavelet approachwhich can be applied to study of diffusions.
Acknowledgments
This work was supported by National Academy of Sciences of Ukraine under ProjectVC-138.
Appendix