Numerical Solution to the Drift-Diffusion Equation of Point

Document Sample
Numerical Solution to the Drift-Diffusion Equation of Point Powered By Docstoc
					Numerical Solution to the Drift-Diffusion Equation of Point Defects under Spherical Symmetry
Mohammed Suleiman Hussein Suleiman ( African Institute for Mathematical Sciences (AIMS)
Supervised by: Prof. D.T. Britton and Prof. M. Harting University of Cape Town, South Africa

22 May 2008
Submitted in partial fulfillment of a postgraduate diploma at AIMS

The diffusion of point defects is considered as a basis for reactions in solids which change important physical properties of materials. In this work, after reviewing some essential concepts on the subject, a first approximation to the problem of clustering of point defects in cubic crystallines has been made by considering certain spherical symmetries. Then, a numerical solution to the obtained drift-diffussion equation is acquired using a finite difference technique. The obtained results are discussed within the conceptual framework introduced.

I, the undersigned, hereby declare that the work contained in this essay is my original work, and that any work done by others or by myself previously has been acknowledged and referenced accordingly.

Mohammed Suleiman Hussein Suleiman, 22 May 2008 i

Abstract 1 Introduction 2 Crystalline Solids 2.1 2.2 2.3 2.4 2.5 2.6 Space Lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Unit and Primitive Cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Crystal Systems and Bravais Lattices . . . . . . . . . . . . . . . . . . . . . . . . . . . . Interstices in Crystal Structures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Crystal Binding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Crystal Defects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . i 1 2 2 2 3 3 4 5 7 7 7 7 9 10

3 Diffusion and Diffusion Equations 3.1 3.2 3.3 3.4 Fick’s First Law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fick’s Second Law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Forced Diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The Fokker-Planck (FP) Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4 Diffusion in Crystalline Solids 4.1 4.2 4.3 4.4

Diffusivity of Crystalline Solids . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 D Tensor of Cubic Crystals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 D Tensor of Non-Cubic Crystals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

Mechanisms of Diffusion in Crystalline Solids . . . . . . . . . . . . . . . . . . . . . . . 12 14

5 Diffusion of Point Defects in Cubic Crystals Under Spherical Symmetry 5.1 5.2 5.3 5.4

Applicability of the FP Equation to the Problem of Point Defects Transport in CubicCrystals . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 Spherical Symmetries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 The FP Equation in Spherical Polar Coordinates . . . . . . . . . . . . . . . . . . . . . . 15 Numerical Solution to the FP Equation . . . . . . . . . . . . . . . . . . . . . . . . . . 16 20

6 Results and Discussion 6.1

Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 ii

6.1.1 6.1.2 6.1.3 6.2

Results for True Diffusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Results for Diffusion with an Outward Driving Force . . . . . . . . . . . . . . . . 21 Results for Diffusion with an Inward Driving Force . . . . . . . . . . . . . . . . . 21

Results Justification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21 23 25

7 Conclusion References


1. Introduction
The subject of diffusion in solids has been an active area of research for a long time. Physical properties of materials are strongly related to their internal structure. A real understanding of changes of many physical properties must be based on knowledge of diffusion, since most changes in the structure of materials occur by diffusion [She63]. The purpose of this work is to introduce some basic notions of diffusion in crystalline solids, then to try to apply these essentials to model and solve a problem of early stage clustering of point defects in a crystalline solid bulk. The treatment is very general, that is, no real system is considered. A basic assumption and some simplifications are made to solve the problem. Thus, no attempt is made at completeness or rigor. Rather it is a first step or approximation. Since crystalline solids are to be considered, in our problem, as the mediums of diffusion, their internal geometry is discussed in chapter 2. The nature of the diffusing particles, point defects, are also included. In chapter 3, the general phenomenon of diffusion of particles through a medium is introduced briefly and in a very general sense. That is, in a way that the general derived differential equation can be equally applied to any physical diffusion phenomenon. Diffusivity D and mobility η are two basic quantities which contain the physics of any diffusion system [BH02]. Their mathematical nature is given in chapter 4. Then this will be followed by a brief discussion about how a particle diffuses in a crystal lattice. Chapter 5 starts with introducing the basic assumption of this work: the diffusion of point defects can be treated by ignoring the crystalline structure of the medium. Then, to simplify the problem, three spherical symmetries are assumed. According to this assumption and to the simplifications made, it is shown that the problem becomes a matter of solving the obtained partial differential equation. Then, because of the geometry of the problem, this obtained equation is rewritten in spherical polar coordinates, and a numerical method of solution is investigated. The obtained results are presented and justified in chapter 6, followed by a conclusion and recommendations for future work in the last chapter.


2. Crystalline Solids
According to their internal geometry, solids can be classified into two main types: crystalline and noncrystalline or amorphous [All01]. Each type is characterized by the size of an ordered region within the material. In crystals, atoms or molecules have a very regular geometric periodicity over extended regions, i.e. over a large number of atoms or molecules. In contrast, amorphous materials have order only over a few atomic or molecular dimensions [Tab00][Nea03]. Structural and chemical defects exist in both structures, but for clarity the discussion below is focused on defects in crystals.


Space Lattice

A crystal is constructed by repetition, in the three dimensions, of identical atoms or groups of atoms which are called bases. Then we define the lattice, or space lattice, which is a mathematical abstraction, as the set of mathematical points to which bases are attached [Kit05].


Unit and Primitive Cell

A small parallelepiped of the crystal that can be used to to produce the entire crystal is known as a unit cell. It is not a unique entity. The smallest unit cell is called the primitive cell. In many cases, it is more convenient to use a unit cell that is not primitive[Nea03].A generalized primitive unit cell is shown in figure 2.1.

Figure 2.1: A generalized primitive unit cell. Modified from: [Nea03]

The relationship between the unit primitive cell and the lattice can be given by the three primitive 2

Section 2.3. Crystal Systems and Bravais Lattices

Page 3


γ α β
b a

Figure 2.2: A generalized unit cell.

translation vectors: a1 , a2 and a3 , such that every two lattice points r and r in the crystal can be related to each other by ′ r = r + n1 a1 + n2 a2 + n3 a3 , (2.1) where n1 , n2 and n3 are integers[Nea03][Kit05]. So, the set of points r given by 2.1 for all n1 , n2 and n3 defines the lattice, and the volume of the unit primitive cell is given by [Kit05]: v = a1 · a2 × a3 . (2.2)



Crystal Systems and Bravais Lattices

Crystal lattices can be carried into themselves by the lattice translation T = n1 a1 + n2 a2 + n3 a3 and by various other symmetry operations. A typical symmetry operation is that of the rotation about an axis that passes through a lattice point. Lattices can be found such that one-, two-, three-, four- and six-fold rotation axes carry the lattice into itself, crresponding to rotation by 2π, 2π/2, 2π/3, 2π/4 and 2π/6 radians and by integral multiples of these rotations. For example, the a cubic lattice has three four-fold rotation axes, four three-fold rotation axes and six two-fold rotation axes [Kit05]. It was found that - by applying equation 2.1 and the other possible symmetry operations - there are only seven ways to draw a unit cell. These possible cells are called crystal systems, and when points are distributed on these systems it produces fourteen possibilities called Bravais lattices. Each cell is characterized by its axial relations, that is the lenghts of its three axes a, b, c and the angles α, β, γ between them as shown in figure 2.2 [All01]. The Bravais unit cells are given in table 2.1 [Kit05]. The three cubic crystals are the simple cubic (sc) lattice, the body-centered cubic (bcc) and the facecentered cubic (fcc) lattice [Kit05]. Their unit cells are shown in figure 2.3 and they will be considered in chapter 4.


Interstices in Crystal Structures

By knowing the crystal structure of a material and its lattice dimensions, several charateristics of the crystal can be determined such as the volume density of atoms. Two factors are always used: (i) the packing density is defined as number of atoms per unit cell [Nea03] (ii) the packing fraction is the

Section 2.5. Crystal Binding Table 2.1: Crystal Systems Restrictions on Axes’s Lengths Restrictions on Angles a=b=c α=β=γ a=b=c α = β = 90◦ = γ a=b=c α = β = γ = 90◦ a=b=c α = β = γ = 90◦ a=b=c α = β = γ = 90◦ a=b=c α = β = γ = 90◦ a=b=c α = β = 90◦ γ = 120◦

Page 4

System Triclinic Monoclinic Orthorhombic Tetragonal Cubic Trigonal Hexagonal

Number of Lattices 1 2 4 2 3 1 1

Figure 2.3: The three cubic cells: (a) simple ubic, (b) body-centered cubic, (c) face-centered cubic. Source: [Nea03]

maximum proportion of the available volume -of the unit cell- that can be filled with hard spheres [Kit05]. The unoccupied spaces in the lattice are called interstices. Interstices have certain geometrical shapes and certain number in a unit cell [All01]. For example, the interstitial sites of an fcc lattice also form an fcc lattice [She63]p.43. The existence of the interstices in the crystal play important roles such as hosting the interstitial atoms in the interstitial solid solutions [All01] and in the diffusion of point defects.


Crystal Binding

The cohesive energy of a crystal is defined as the energy required to seperate the crystal components into neutral free atoms at rest, at infinite seperation, with the same electronic configuration [Kit05].The types of bonds in crystals are: 1. Ionic bonds: Oppositely charged ions experience a coloumb attraction and form a bond refered to as ionic bond. If the ions were to get too close, a repulsive force would become dominant, so an equilibrium distance results between these two ions [Nea03]. The term lattice energy is used in the discussion of ionic crystals and is defined as the energy that must be added to the crystal to seperate its component ions into free ions at rest at infinite seperation [Kit05].

Section 2.6. Crystal Defects

Page 5

2. Covalent bonds: Here atoms tend to achieve closed-valence energy shells by sharing their outer electrons [Nea03][FLS75]. 3. Metallic bonds: Positive metallic ions might be thought of as being sorrounded by a sea of negative electrons, the solid being held together by the electrostatic forces [Nea03][FLS75]. 4. Van der waals bonds: In some molecules the effective centre of the positive charge of the molecule is not the same as the effective centre of the negative charge. This nonsymmetry in the charge distribution results in a small electric dipole that can interact with a dipole of other molecules [Nea03]. In many materials the bonds between atoms in a crystal are mixed. That is, bonds are not purely ionic, purely covalent or purely metallic. Instead, it is a mixture of two or more bond types [All01]. The different strength of some crystals along different directions depends on the kind of interatomic binding in the different directions [FLS75].


Crystal Defects

The term defect can be interpreted to mean any equilibrium deviation from the perfect homogeneity and the geometric periodicity of the crystal structure. According to this definition, fluctuations are not considered. Crystal defects can be classified, according to their geometrical shape and according to the dimensions that they take, into four main classes [All01]: 1. Point Defects: The word point is used to denote the localized nature of the defects. Theoretically, point defects are zero dimensional, but in reality involve single atoms or single-atomic locations. The commom point defects in crystals are: (i) foreign atoms, or impurities, may be located at normal lattice sites or between normal sites and they are called substitutional impurities and interstitial impurities, respectively [Nea03], (ii) vacancies, that is vacant lattice sites and (iii) self-interstitials, which are extra atoms of the same type of the crystal but not in regular lattice positions. Although point defects are generally introduced by non-equilbrium processes, such as ion implantation [HNR+ 07] or irradiation [PHI91], there is always thermal equilibrium concentrations of point defects present. If the thermal energy of an atom in a crystal exceeds the energy required to remove an atom from its lattice site then that atom can leave its site. As a result, a number of point defects can be formed in the structure: (a) Isolated vacancy: If the atom released and migrates to the surface, or to a grain boundary, it will leave an isolated vacancy inside the crystal. This defect is also known as a Schottky defect, and is illustrated in figure 2.4(a). The common Schottky defect occurs in ionic crystals, where two vacancies are formed in two different lattice sites because of the imigration of a positive ion “cation” and a negative ion “anion”. This type of defect is also called the ion-pair defect, in which the crystal is roughly kept electrostaticaly neutral on a local scale. (b) Isolated interstitial: An atom can diffuse in from the surface of the crystal and transfer to an interstitial position, that is a position not normally occupied by an atom, forming an isolated interstitial (figure 2.4(b)).

Section 2.6. Crystal Defects

Page 6




Figure 2.4: Two dimensional illustration of the formation of: (a) an isolated vacancy, (b) an isolated interstitial, and (c) a vacancy-interstitial pair. (c) Frenkel defect: If the atom released and migrates only a short distance from the vanacy that it created, the two form an interstitial-vacancy pair called a Frenkel defect, as illustrated in figure 2.4(c). In higly-packed structures the energy required to introduce an interstitial atom of the element is so large and that make their occurence extremely unlikely. In contrast to that, interstitials are more easily accommodated in structures with low packing fracture [Wal92][All01][Kit05]. 2. Line Defects (Dislocations): The point defects involve single atoms or single-atoms locations. In contrast, dislocation is a deviation from the perfect lattice that involve an entire row of atoms or atoms sites. There are many types of dislocations such as edge and screw dislocation. Dislocations affect most of the physical properties of the crystal [Nea03] [All01]. 3. Planar Defects: The previously described crystals are known as single-crystal materials, that is, they have a high degree of order and regular geometric periodicity throughout the entire volume of the material. In other materials the regular geometric periodicity may extend only over many atomic or molecular dimensions forming single-crystal regions which usually vary in size and in orientation with respect to one another. These single-crystal regions are called grains and they are seperated from one another by surfaces called grain boundaries [Nea03][All01]. 4. Bulk Defects: When vacancies cluster together they form “voids”, and when foreing atoms cluster together they are called “precipitates”. These clusters are called bulk defects. That is, bulk defects are clusters of point defects. These defects affect the physical properties of the perfect crystals, and they can interact with each other forming complex defects [Kit54][Nea03].

3. Diffusion and Diffusion Equations
Although later the treatment will be restricted specifically to the diffusion of certain speicies in certain mediums, in this chapter, the general phenomenon of diffusion of particles through a medium will be briefly introduced. That is, the following discusion can be equally applied to gas diffusion and charge transport in semiconductors. Hence, the term particle will be used in a very general sense, which includes even crystal vacancies.


Fick’s First Law

True diffusion, or simply diffusion, of a certain species in a medium is the outcome of a random motion which, by virtue of a concentration gradient of that species results in a net flux of that species away from the highest concentration. The net diffusive flux JD is related to the gradient of the concentration C(x, t) by a famous phenomenological relation called Fick’s first law: JD = −D∇C(x, t) , (3.1)

where C(x, t) is the concentration of particles at a position x and at time t, (∇C) is the concentration gradient, and D is called the diffusion coefficient or diffusivity, which is defined as the net flux of particles per unit concentration gradient [She63] [BH02]. If there are other types of gradient in the system, other terms are added to Fick’s flux equation [She63] [BH02]. This will be discussed later.


Fick’s Second Law

Equation 3.1 can only be used in the steady state of flow, that is, when the net flux J is independent of time: J = −D∇C(x, t) = const. (3.2) In the time-dependent case, the rate of change of particles concentration at any point x in the distribution can be related to the net flux J by ∂ C(x, t) = −∇ · J . (3.3) ∂t Equation 3.3 is known as Fick’s second law and can be derived from equation 3.1. It simply states that a net flow out of a volume results in a decreased density inside the volume. This equation appears in a wide variety of physical problems, ranging from a probability density in quantum mechanics to neutrons leakage in a nuclear reactor [Arf85] [She63].


Forced Diffusion

The flow of particles can in fact be due to two causes: one is the mentioned concentration gradient, and the other is the action of a driving force. That is, if we consider a single particle moving in a potential


Section 3.3. Forced Diffusion field V (x); the potential gradient exerts a force F on the particle given by the equation F = −∇V (x) .

Page 8


It is found empirically that under the influence of this force, particles move with a certain average terminal elocity vd . This fact is mathematically expressed by the equation vd = ηF , (3.5)

where η, is called mobility, and can be defined as the ratio between the average terminal velocity of the particles and the applied force. It is a measure of the response of the particles of the medium to an applied force [She63] [PHI91]. For a concentration distribution of particles C(x, t) under the influence of a driving force F(x), vd gives rise to a flux: JF = C(x, t)vd = C(x, t)ηF. (3.6) Accordingly, one can define the forced diffusion as a directed motion that is superimposed onto the random motion of the particles. Forced diffusion is also refered to as drift [BH02]. Three kinds of driving force that can drive particles, including point defects, in a solid medium are briefly mentioned below. Each of them can be exspressed as a potential gradient [Wip76][All01] [PHI91]: 1. Electric field: The motion of charged particles under an electric field E is called electroimigration or electrotransport. The electric force Fe is given by Fe = q ∗ E , (3.7)

where q ∗ is the effective electric charge of the particle, and E can be expressed in terms of the electrical potential Φ as E = −∇Φ(x) . (3.8) 2. Temperature gradient: A temperature gradient will result in a force given by Fth = − Q∗ ∇T (x) , T (3.9)

where Q∗ is called migration enthalpy or heat of transport, and the motion of particles under temperature gradient is called thermoimigration or thermotransport. 3. Stress gradient: Stress gradient will result in a force given by Fσ = −∇U (x) , (3.10)

where U is the elastic interaction energy in the stress field σ(x). A uniform stress cannot generate a flux (Curie-Neumann principle); thus the driving force must be the stress gradient and its effect must be considered whenever the elastic interaction energy of the particle with the stress field is large enough. A typical case is Gorski effect: imigration of atoms in an elastically deformed sample.

Section 3.4. The Fokker-Planck (FP) Equation

Page 9


The Fokker-Planck (FP) Equation

The total net flux J in equation 3.3 is due to both contributions of JD and JF defined by equations 3.1 and 3.6 respectively. That is J = JD + JF = −D∇C(x, t) + ηC(x, t)F. Substiuting for J from equation 3.11 into equation 3.3 yields ∂ C(x, t) = ∇ · (D∇C(x, t)) − ∇ · (ηC(x, t)F). ∂t (3.12) (3.11)

Equation 3.3, or equation 3.12, is a continuity equation, that is, it is valid only for species which obey a conservation law [PHI91]. If one deals with entities which are not conserved, then an additional term S(x, t) which is equal to the rate of production or destruction of these entities per unit volume should be added to the right hand side of equation 3.3, or equation 3.12 [She63] [BH02]. Thus ∂ C(x, t) = −∇ · J + S(x, t) , ∂t or (3.13)

∂ C(x, t) = ∇ · (D∇C(x, t)) − ∇ · (ηC(x, t)F) + S(x, t) . (3.14) ∂t Equation 3.14 is the general form of the drift-diffusion equation [BH02], and it is known in the literature as a Fokker-Planck (FP) equation1 [GT88].

1 A Fokker-Planck (FP) equation is a partial differential equation for the time evolution of the probability density of a particle (or field). That is, for a diffusing particle, one writes a differential equation for C(x, t), where C(x, t) is to be defined as the probability that the diffusing particle is to be found in the interval x + dx at time t [LMN07]. Equation 3.14 is one of the simplest Fokker-Planck equations [Ris96].

4. Diffusion in Crystalline Solids
While in chapter 3 the diffisivity D was defined and its units were mentioned, but nothing more about the mathematical nature of D or its dependence on the structure of the meduim was mentioned. This will be given in the first sections of this chapter in order to gustify the applicability of the diffusion equations introduced in chapter 3 to the problem of the diffusion in crystalline solids. Then this will be followed by a brief discussion about how a particle diffuses in a lattice.


Diffusivity of Crystalline Solids

In equation 3.1, the two vectors J and ∇C are, in the most general case, not parallel, and then D cannot be considered as a constant. Instead, for the general case, it is assumed that a given component of the flux is influenced by each of the components of the gradient. Thus, in the cartesian coordinates, this situation is described by the equations: Jx = −D11 Jy Jz ∂C ∂C ∂C − D12 − D13 ∂x ∂y ∂z ∂C ∂C ∂C − D22 − D23 = −D21 ∂x ∂y ∂z ∂C ∂C ∂C − D32 − D33 = −D31 dx dy ∂z , , . (4.1) (4.2) (4.3)

where the three scalar fluxes Jx , Jy and Jz are parallel to the three cartesian coordinates, and the set of nine numbers designated Dij is called a second-order tensor and is defined by the above equations. The same argument can be made for the mobility η in equation 3.6 [She63][FLS75].

The above system of equations can be rewritten in the matrix form as      D11 D12 D13 ∂C/∂x Jx −  D21 D22 D23   ∂C/∂y  =  Jy  D31 D32 D33 ∂C/∂z Jz




D Tensor of Cubic Crystals

Let us consider a cube of a material whose lattice is a cubic one. Taking a single crystal and setting up its axes parallel to the cartesian coordinates, and a concentration gradient is established such that ∂C/∂x = α = 0 and ∂C/∂y = ∂C/∂z = 0. Equation 4.1 then gives the component of the flux parallel to the gradient as ∂C . (4.5) Jx = −D11 ∂x If this gradient is removed and replaced by an equal gradient along the y axis, that is, ∂C/∂y = α = 0 and ∂C/∂x = ∂C/∂z = 0, then the component of the flux parallel to the gradient will be Jy = −D22 ∂C ∂y . (4.6)


Section 4.2. D Tensor of Cubic Crystals The same argument can be made for Jz and Jz = −D33 ∂C ∂z .

Page 11


In a cubic lattice crystal x, y and z axes are indistinguishable. And since the gradient were of equal magnitudes in the three cases, the cubic symmetry requires that Jx = Jy = Jz in the three experiments. This implies that D11 α = D22 α = D33 α, that is D11 = D22 = D33 . (4.8)

Again, let us establish a gradient such that dC/dy = α = 0 but ∂C/∂x = ∂C/∂z = 0. By measuring the component of the flux in the x axis, one can determine D12 , since equation 4.1 gives Jx = −D12 ∂C = −D12 α . ∂y (4.9)

Now, let us imagine that the source and sink which established this gradient are removed and the crystal rotated 180◦ about its x axis, relative to the source and sink. Since the lattice x direction has four-fold rotational symmetry, the original and final positions of the lattice will be indistinguishable. If the source and the sink are again applied, the gradient is again established along the y axes. However, since the y axis is fixed in the crystal, the 180◦ has interchanged the +y and the −y axes, so the gradient will be negative of its previous value. This means that now ∂C/∂y = −α, and ∂C/∂x = ∂C/∂z = 0. This gives ∂C = D12 α . (4.10) Jx = −D12 ∂y But by symmetry, Jx from equation 4.9 must equal Jx from equation 4.11, that is D12 α = −D12 α . (4.11)

or more simply

This requirement of symmetry can be satisfied only if D12 = 0. If the gradient had been along the z axis, the same rotation could have been used to show that D13 = 0. Since the y and the z axes also have two-fold rotation symmetry, it follows that all off-diagonal terms are zero. Accordingly, for cubic crystals   D11 0 0 D =  0 D11 0  , (4.12) 0 0 D11 D = D(x). (4.13)

The same argument can be made for the mobility η in equation 3.6, that is η = η(x). (4.14)

for cubic crystals [She63] [FLS75]. Equations 4.13 and 4.14 simply show that cubic crystals are isotropic.

Section 4.3. D Tensor of Non-Cubic Crystals

Page 12


D Tensor of Non-Cubic Crystals

and for orthohrombic crystals

It is shown in reference [She63] that equation 4.12 is no longer true for non-cubic crystals. For example, for tetragonal and hexagonal lattices D is given by   D11 0 0 D =  0 D11 0  , (4.15) 0 0 D33  D11 0 0 D =  0 D22 0  0 0 D33  . (4.16)

That is, non-cubic crystals are inisotropic.


Mechanisms of Diffusion in Crystalline Solids

In both cases, true and forced diffusion, for a particle to diffuse in a crystal, it must overcome the potential energy barrier presented by its neighbours. If the barrier height is E and the particle characteristic vibrational frequency is ν, then the propability ρ that during a unit time this particle will have enough thermal energy to pass over the barrier is ρ = νe−E/kT , (4.17)

that is, each unit time the particle makes ν jumps, with a probability e−E/kT of overcoming the barrier in each try. ρ is called the jump frequency [Kit05]. Most mechanisms of diffusion require the presence of point defects discussed in section 2.6. There are three possible elementary mechanism of diffusion of particles in a solid meduim [Kit05][Wal92]: 1. Ring Mechanism: For an atom in a lattice site, it is easy to conceive a mechanism of atomic motion by direct interchange by a rotation about a midpoint (figure 4.1a), or by rotation of a ring of three or more atoms 1 [PHI91][She63]. 2. Interstitial Mechanism: An atom is jumping from interstitial site (discussed in section 2.4) to interstitial site is a point defect, and is said to diffuse by interstitial mechanism. This is illustrated, in two dimensional lattice, in figure 4.1b. An appreciable local dilation of the lattice must occur before the jump can occur. It is this dilation or distortion which constitutes the barrier to an interstitial atom changing sites. A variety of interstitial mechanisms is possible. In the direct interstitial mechanism, the atom passes without permanently displacing any of the matrix atoms and also there is no dependence between the successive jumps. In contrast, in the indirect interstitial mechanism (also known as interstitialcy mechanism) a relatively large interstitial atom pushes one of its nearest-neighbour atoms into an interstitial position and occupied the lattice site previously occupied by the displaced atom [PHI91] [She63].
The ring mechanism is not known to operate in any metal or alloy, but it has been suggested as a mechanism which would explain some apparent anomalies in D for bcc metals [She63].

Section 4.4. Mechanisms of Diffusion in Crystalline Solids

Page 13

3. Vacancy Mechanism: If a lattice site is not occupied, a nearest neighbour atom can jump onto this site, and the vacancy will appear on the site that the atom has just vacated, the atom is said to have diffused by a vacancy mechanism [PHI91][She63]. This mechanism is illustrated in figure 4.1c in two dimensional lattice.


( b)

( c )

Figure 4.1: Three basic mechanisms of diffusion : (a) Interchange by rotation about a common midpoint. (b) Migration through interstitial sites. (c) Exchange of position with a vacancy.

More complex mechanisms are possible, and several elementary mechanisms can be operative simultaneously. And if different types of particles participate in the diffusion process (e.g. host atoms, vacancies and interstitials), then the succesive jumps of the particles are dependent [PHI91].

5. Diffusion of Point Defects in Cubic Crystals Under Spherical Symmetry
In this chapter, the problem of diffussion of point defects in cubic crystals will be treated. After introducing some basic assumptions to simplify the problem, the Fokker-Planck equation, equation 3.14, will be rewritten, for such homogenous isotropic solid mediums, in spherical coordinates. Then a numerical method of solution, given an initial distribution C(x, 0) and defined driving fields F(x), will be investigated.


Applicability of the FP Equation to the Problem of Point Defects Transport in Cubic-Crystals

Fick’s law (equation 3.1), continuity equation (equation 3.3) and the consequent derived equations are valid for continuum mediums and continous C(x, t) [PHI91]. Fick treated the diffusion problem by assuming that the medium is continuum, neglecting the discontinuity of the real material structure [PHI91] [She63] and by ignoring the diffusion mechanisms [She63]. We are going to make the same assumption since all subsequent developments in the theory of solids have in no way affected the validity of Fick’s approach [She63]. Accordingly, under this assumption, at no point in the following discussion the structural nature of the diffusion medium or the diffusion mechanism will enter the problem. Moreover, for cubic crystals, equations 4.13 and 4.14 become D = const. , η = const. That is, neither D nor η has spatial dependence, and the medium is said to be homogeneous. (5.1) (5.2)


Spherical Symmetries

To simplify the problem more, we consider a single cubic crystal with a spherical boundary, or a polycrystalline material which has cubic structure, then imagining one of the grains has a spherical boundary. Also a Gaussian initial distribution is assumed, that is C0 = C(r, θ, φ, 0) = N e− for some constant N 1 , and with boundary condition: C(R, θ, φ, 0) = N e− C(R, θ, φ, t) = 0,
(R−r0 )2 2σ 2 (r−r0 )2 2σ 2




t=0 t > 0,

(5.4) (5.5)

where R is the radius of the bulk or the grain, 0 < r0 < R is the expected value with 0 < r0 < R, and σ is the standard deviation.
The prefactor N can be chosen as a normalization factor in order for C(r, θ, φ, 0) to satisfy the definition of the FP equation (equation 3.14). That is C(r, θ, φ, t) must be redefined as a probability density.


Section 5.3. The FP Equation in Spherical Polar Coordinates

Page 15

Also, the case when there is only true diffusion will be considered first, then a driving force with spherical symmetry, that is 1 F ∼ ± 2 er (5.6) r will be introduced.


The FP Equation in Spherical Polar Coordinates

Equation 3.14 yields ∂ C = D∇2 C + ∇D · ∇C − C∇η · F − η∇C · F − Cη∇ · F + S ∂t . (5.7)

And since we assumed that the medium which we are dealing with is homogenous and isotropic, then both ∇D and ∇η will vanish, since D and η are constants, and equation 5.7 becomes ∂ C = D∇2 C − η(∇C) · F − Cη∇ · F + S ∂t . (5.8)

Also, since the bulk or the grain boundary, the driving force and the initial distribution all have spherical symmetry, it is better then to write equation 5.8 in spherical polar coordinates. In this coordinates we have [Arf85]: ∇C(r, θ, φ, t) = er ∇2 C(r, θ, φ, t) = + ∇ · F(r, θ, φ) = 1 ∂ 1 ∂ ∂ C(r, θ, φ, t) + eθ C(r, θ, φ, t) + eφ C(r, θ, φ, t) , ∂r r ∂θ rsinθ ∂φ ∂ ∂ ∂ 1 ∂ (sinθ C(r, θ, φ, t)) (5.9) sinθ (r 2 C(r, θ, φ, t)) + 2 sinθ r ∂r ∂r ∂θ ∂θ 1 ∂2 C(r, θ, φ, t) , sinθ ∂φ2 ∂ ∂ ∂ 1 . sinθ (r 2 Fr (r, θ, φ)) + r (sinθFθ (r, θ, φ)) + r Fφ (r, θ, φ) r 2 sinθ ∂r ∂θ ∂φ

Substituting for ∇C(r, θ, φ, t), ∇2 C(r, θ, φ, t) and ∇ · F(r, θ, φ) from 5.9 into equation 5.8 yields

∂ C(r, θ, φ, t) = D ∂t

∂2C ∂2C 2 ∂C 1 ∂2C cotθ ∂C 1 + + 2 2 + 2 + 2 2 ∂r 2 r ∂r r ∂θ r ∂θ r sin θ ∂φ2 1 ∂C 1 ∂C ∂C + eθ + eφ · (er Fr + eθ Fθ + eφ Fφ ) − η er ∂r r ∂θ rsinθ ∂φ ∂Fr 2 1 ∂Fθ cotθ 1 ∂Fφ − Cη + Fr + + Fθ + ∂r r r ∂θ r rsinθ ∂φ + S(r, θ, φ, t) .


This is the FP equation in the spherical coordinates, which must be solved for a given initial distribution and a well realistic defined driving force F to model the clustering of point defects.

Section 5.4. Numerical Solution to the FP Equation

Page 16


Numerical Solution to the FP Equation

Except for a few special cases, there is no straightforward analytical solutions to such drift-diffusion equation [Bri94] [Ris96]. Generally, however, it is difficult to obtain analytical solutions especially if no seperation of variables is possible, and its numerical solution is a non-trivial problem in numerical analysis [GT88]. The clearest numerical method is to follow the dynamical processes step by step using a finite differencing technique [Bri94]. In this technique, first we represent C(r, θ, φ, t) by its values at discrete set of points: ri = i△r, θj = j△θ, φk = k△φ, tn = n△t, i = 0, 1, 2, ..., I j = 0, 1, 2, ..., J k = 0, 1, 2, ..., K n = 0, 1, 2, ..., N (5.11)

where △r, △θ and △φ are the chosen spatial spacings, and △t is the temporal one[Pea02].

The Explicit (FTCS) Method:
The simplest finite difference method is to calculate C(r, θ, φ, t) explicitly. That is, at each time step C(r, θ, φ, t) can be calculated explicitly from its values in the previous time step [Pea02]. In this scheme, the numerical form of equation 5.10 can be written as

Section 5.4. Numerical Solution to the FP Equation

Page 17

C(r, θ, φ, t + △t) − C(r, θ, φ, t) △t

C(r + △r, θ, φ, t) − 2C(r, θ, φ, t) + C(r − △r, θ, φ, t) (△r)2 2 C(r + △r, θ, φ, t) − C(r − △r, θ, φ, t) + r 2△r 1 C(r, θ + △θ, φ, t) − 2C(r, θ, φ, t) + C(r, θ − △θ, φ, t) + r2 (△θ)2 cotθ C(r, θ + △θ, φ, t) − C(r, θ − △θ, φ, t) + r2 2△θ C(r, θ, φ + △φ, t) − 2C(r, θ, φ, t) + C(r, θ, φ − △φ, t) 1 + 2 sin2 θ r (△φ)2 C(r + △r, θ, φ, t) − C(r − △r, θ, φ, t) −η Fr (r, θ, φ) 2△r Fθ (r, θ, φ) C(r, θ + △θ, φ, t) − C(r, θ − △θ, φ, t) + r 2△θ Fφ (r, θ, φ) C(r, θ, φ + △φ, t) − C(r, θ, φ − △φ, t) + rsinθ 2△φ Fr (r + △r, θ, φ) − Fr (r − △r, θ, φ) 2 −η C(r, θ, φ, t) + Fr (r, θ, φ) 2△r r Fθ (r, θ + △θ, φ) − Fθ (r, θ − △θ, φ) cotθ 1 + Fθ (r, θ, φ) + r 2△θ r 1 Fφ (r, θ, φ + △φ) − Fφ (r, θ, φ − △φ) + rsinθ 2△φ + S(r, θ, φ, t). (5.12) = D

n If we write Cijk for C(ri , θj , φk , tn ), Fr,ijk for Fr (ri , θj , φk ) ... and so on, equation 5.12, in this notations, becomes

Section 5.4. Numerical Solution to the FP Equation

Page 18

n+1 n Cijk − Cijk


= + + −η +


n n n C(i+1)jk − 2Cijk + C(i−1)jk


n n 2 C(i+1)jk − C(i−1)jk + ri 2△r

n n n n n cotθj Ci(j+1)k − Ci(j−1)k 1 Ci(j+1)k − 2Cijk + Ci(j−1)k + 2 2 (△θ)2 2△θ ri ri n n + Cn Cij(k+1) − 2Cijk 1 ij(k−1) 2 ri sin2 θ

(△φ)2 2△r +
n n Fθ,ijk Ci(j+1)k − Ci(j−1)k ri 2△θ


n n C(i+1)jk − C(i−1)jk

n n Fφ,ijk Cij(k+1) − Cij(k−1) ri sinθj 2△φ

n −η Cijk

+ +

Fr,(i+1)jk − Fr,(i−1)jk 2Fr,ijk + 2△r ri cotθj Fθ,ijk 1 Fθ,i(j+1)k − Fθ,i(j−1)k + ri 2△θ ri 1 Fφ,ij(k+1) − Fφ,ij(k−1) n + Sijk ri sinθj 2△φ


This finite difference approximation to equation 5.10 is known as FTCS representation (Forward Time Centered Space). It is an easy algorithm to derive, takes little storage in the computer and executes quickly [Pea02]. But this method is said to be unstable if round off errors or other errors grow rapidly as the computations proceed. The necessity of using very small time steps is the principal fault of this method [ZC97]. Following this scheme, a code has been written in Python programming language to solve equation 5.13 n+1 for Cijk at each time step. The concentration at each spatial point is calculated as follows: 1. Using 5.11 the space and time were configurated with N = 100, I = 60, J = 50 and K = 50. 2. An N × I × J × K matrix was created to store the computed values of C 2 . 3. The radius R of the bulk and the total time T of the observation were both chosen to be unity for some arbitrary units. 4. r = 0 was avoided, since it represents a singularity in equation 5.10. 5. The values D = 1 × 10−2 and η = 1 × 10−4 were chosen. 6. The initial distribution 5.3 was stored with an expected value at r0 = 0.4R and standard deviation σ = 0.1R.
n 7. The Dirichlet boundary condition: CIjk = 0 for n > 0 was introduced. n 8. Neumann boundary conditions were applied for the angular dependence. That implies: Ci0k = n = C n , and C n n n , Cn = Cn Ci1k iJk ij1 ijK = CijK−1 . i(J−1)k , Cij0

9. Radial concentration distributions at some selected time steps were plotted.

It could be thought of as N matrices each of order I × J × K.

Section 5.4. Numerical Solution to the FP Equation

Page 19

The Implicit (BTCS) Method:
In this scheme the the spatial derivatives on the right-hand side are evaluated at time step n + 1. That is
n+1 n Cijk − Cijk


= + + −η +


n+1 n+1 n+1 C(i+1)jk − 2Cijk + C(i−1)jk


n+1 n+1 2 C(i+1)jk − C(i−1)jk + ri 2△r

n+1 n+1 n+1 n+1 n+1 cotθj Ci(j+1)k − Ci(j−1)k 1 Ci(j+1)k − 2Cijk + Ci(j−1)k + 2 2 (△θ)2 2△θ ri ri n+1 n+1 n+1 Cij(k+1) − 2Cijk + Cij(k−1) 1 2 (△φ)2 ri sin2 θ


n+1 n+1 C(i+1)jk − C(i−1)jk



n+1 n+1 Fθ,ijk Ci(j+1)k − Ci(j−1)k ri 2△θ

n+1 n+1 Fφ,ijk Cij(k+1) − Cij(k−1) ri sinθj 2△φ

n+1 −η Cijk

+ +

Fr,(i+1)jk − Fr,(i−1)jk 2Fr,ijk + 2△r ri Fθ,i(j+1)k − Fθ,i(j−1)k cotθj Fθ,ijk 1 + ri 2△θ ri 1 Fφ,ij(k+1) − Fφ,ij(k−1) n+1 + Sijk ri sinθj 2△φ


Such a scheme is called fully implicit or backward time. To solve equation 5.14 one has to solve a set n+1 of simultaneous linear equations at each time step for Cijk [Pea02]. This method does not suffer from stability problems [ZC97], but, compared to the explicit scheme, it takes more storage in the computer and more time to execute.

6. Results and Discussion
In this chapter, the obtained results are presented, followed by a brief discussion.


Numerical Results

Graphic results for the three cases: F = 0, F = + r12 er and F = − r12 er are shown below. For each obtained figure the computer runs for about 40 minutes.


Results for True Diffusion

Figure 6.1 represents some results obtained in the absence of a driving force. The four different curves show the radial distribution of point defects at different times: (i) at t = 0 the curve represents the given initial Gaussian distribution, (ii) at a later time t = 0.10T the peak’s decrease shows that the diffusion starts to take place, (iii) as the time runs, at t = 0.60T and at t = 0.99T , the decrease of the peak and the broadening of the distribution becomes clearer.

Figure 6.1: Radial distribution of point defects under true diffusion at four different times.


Section 6.2. Results Justification

Page 21


Results for Diffusion with an Outward Driving Force

Figure 6.2 represents some results obtained in the presence of an outward driving force. In addition to the same behaviour shown in figure 6.1 one can notice that: (i) the values of C at r = 0 are greater than those in 6.1 at the same time steps, (ii) there is a formation of a new peak at r = 0.

Figure 6.2: Radial distribution of point defects driven by an outward radial force at four different times.


Results for Diffusion with an Inward Driving Force

Figure 6.3 represents some results obtained in the presence of an inward driving force. In addition to the same behaviour shown in figure 6.1 one can notice that: (i) the values of C at r = 0 are lower than those in figure 6.1 at the same time steps, (ii) there is a formation of a new peak at r = 0.


Results Justification

In all three cases: (i) The decreasing of the concentration at r = r0 implies that some point defects have transported to somewhere else. (ii) The movement of the peak toward the centre r = 0 and the increase of C values at r = 0 show the stage of accumulation of point defects at the centre of our spherical bulk. (iii) The increase in C values in the region r0 << r < R shows that there is also a diffusion in the outward direction but still less than the inward one. (iv) The zero values of C at r = R (at t > 0) satisfies the boundary condition 5.5 and guarantee that the boundary works as a sink for the

Section 6.2. Results Justification

Page 22

Figure 6.3: Radial distribution of point defects driven by an inward radial force at four different times.

diffusing point defects. At r = 0 : (i) There is no concentration gradient at any time, and in the very early stage r = 0 is a stable point. Once the peak at r = 0 is formulated, then it is no longer stable, that is, each particle that reaches r = 0 will be reflected (in our macro scale) by the negative gradient. At this stage, the higher positive gradient in the left side of the old peak will gain and drives more particles toward the centre than the reflected ones. This gives rise to the accumulation at that point. (ii) As the time runs, we expect this peak to keep increasing, while the other two moving peaks1 are decreasing. This will continue untill these two moving peaks diminish and the peak at r = 0 reaches its maximum2 . Then we expect the true diffusion to be totally in the outward direction, enhancing the outward force and retarding the inward one. In the presence of the driving forces: (i) One can notice that the true diffusion behaviour is always dominant. That is, there is no big difference between the velocity of the corresponding peaks in the presence and in the absence of a driving force. This indeed depends on the relative values of D and η, and also on |F|. (ii) We expect the true diffusion toward the centre to be enhanced by the inward force and competed by the outward one. Thus the two behaviours shown in figures 6.2 and 6.3 are the exact opposite of what is expected. It should be because of a minor algebraic error in the code. Unfortunatley, time lacked to check this before the deadline.
1 It is meant by the two moving peaks: the moving one shown in the given figures and the other one on the other half which is not shown here and which should be moving in the opposite direction to the first one. 2 It was not possible to achieve this moment because of the breakdown of the stability of our numerical explicit scheme. That is, the first phase of the dynamics is only studied.

7. Conclusion
In this work: 1. Some essentials on the subject of diffusion in solids have been reviewed. 2. The problem of clustering of point defects has been simplified and modeled succesfully. 3. The algorithm and the achieved code work and give stable numerical solutions within the range of the chosen values for the different variables. 4. In terms of the obtained results, within the conceptual framework introduced in the earlier chapters, and based on the simplifications made in chapter 5: a great partial success has been achieved. 5. Simplifications introduced in chapter 5 allowed to solve the problem and help to gain a better understanding of the problem of clustering of point defects in solids. Thus, the obtained results are to be considered as first approximation.

Future Work
In future work: 1. Asymmetric problems should be considered. That is, asymmetric forces, arbitrary initial distributions and arbitrary boundaries. 2. Realistic forces, diffusivity and mobility values (and forms) should be inserted into the numerical form of the drif-diffusion equation. 3. Non-isotropic (non-cubic) and inhomogenous crystals are to be studied. 4. The stability of the followed explicit scheme must by studied to determine the stability condition for the solution to the FP equation. 5. It is better to follow a stable implicit numerical scheme, which provides more flexibility in applying the realistic values and realistic forms.


“WHO DOES NOT THANK PEOPLE WILL NEVER THANK ALLAH” First of all, I thank my god for the uncountable granted blessings and for the help in completion of this work. Without Him, I would not have finished neither AIMS diploma nor this essay. I am indeed indebted to the guidance I received from my supervisors, Prof. David T. Britton and Prof. Margit Harting, who very patiently provided me with adequate materials and encouraged me to go a head with this work in spite of many set backs. I am grateful to the PhD students of the Solid State and Materials Physics Research Group at UCT: Abdulrafiu Raji, Shadrack Nsengiyumva, and Lloyd Corker, for their assistance and encouragement. I am at loss when I think of expressing my sincere gratitude to the super tutor Matteo Smerlak (a very special man to me) who helped me in this work and during the whole year with his assistance in various forms at various times. To him I owe many things. I also appreciate being helped, in the numerical part of this work, by Dr. Aziz Ouhinou. His help was priceless. My thanks also to everyone gave me a hand with this essay, specially my mates Sara, Felix, Bewketu and Esra. I express my profound gratitude to all AIMS Council, donars and staff; in particular Prof. Neil Turok and Prof. Fritz Hahne for the opportunity that they have given me to study at AIMS. Many special thanks to our father Igsaan Kamalie, Ms. Lynne Teixeira, sister Aida Mpofu, Mr. Jan Groenewald, Ms. Gudrun Schirge, Mr. Emmanuel Kongolo and all Mothers in the kitchen for the wonderful generosity and kindness they showed during our stay at AIMS. I will be failing if I do not mention the angels of AIMS: the lecturers who volunteered themselves for teaching us. Their knowledge, their help, their kindness, their smiles and their patience have left lasting fervent impressions on me. I extend my warm thanks to all our tutors. Special thanks to Dr. Laure Gouba who was always there whenever I needed her. My many thanks are also to our English teacher, Ms. Frances Aron, for her unlimited support. I would like to express my sincere appreciation to my colleagues at AIMS. You guys are all really wonderful. My thanks to all of you for the wonderful atmosphere we shared. May our friendship continue to grow. I wish you all every success in the future. May I ask each and every one of you for forgiveness. And may my god guide all of us to His heaven. I heartily thank my family, to whom I owe everything, for their love, their prayers, their support, their constant encouragement, and for their patience. Finally, I would like to dedicate this essay to my lovely mother and to my two lovely daughters.


[All01] [Arf85] [BH02] [Bri94] [FLS75] [GT88] M. Talb Allah, The science of engineering materials, in Arabic, Omdurman Islamic University, 2001. George Arfken, Mathematical methods for physicsists, third ed., Academic Press, Inc., 1985. D.T. Britton and M. Harting, The influence of strain on point defect dynamics, Advanced Engineering Materials 4 (2002), 629–635. D. T. Britton, Time-dependent diffusion models of the annihilation of low energy positrons implanted in solids, The Royal Society 445 (1994), 57–67. Richard P. Feynman, Robert B. Lighton, and Matthew Sands, The feynman lectures on physics, Addison-Wesley Publishing Company, 1975. Harvey Gould and Jan Tobochnik, An introduction to computer simulation methods, application to physical systems, Addison-Wesley Publishing Company, 1988.

[HNR+ 07] M. Hartin, S. Nsengiyumva, A. T. Raji, G. Dollinger, P. Sperr, S. R. Naido, T. E. Derry snd C. M. Comrie, and D. T. Britton, Near surface stress determination in kr-implanted polycrystalline titanium by x-ray sin2 φ − method, Surface and Coatings Technology 201 (2007), 8237–8241. [Kit54] [Kit05] [LMN07] [Nea03] [Pea02] [PHI91] [Ris96] [She63] [Tab00] [Wal92] [Wip76] [ZC97] Charles Kittel, Introduction to solid state physics, John Wiley and Sons, Inc., 1954. C. Kittel, Introduction to solid state physics, eigth ed., John Wiley and Sons, Inc., 2005. T. B. Liverpool and K. K. Muller-Nedebock, Random walks, polymer and biopolymers: statics and dynamics, second ed., AIMS course study guide, 2007. Donald A. Neamen, Semiconductors physics and devices: Basic principles, third ed., McGrawHill Higher Education, 2003. William H. Press and et al, Numerical recipes in c, the art of scientific computing, 2 ed., Cambridge University Press, 2002. ´ Jean PHILIBERT, Atom movements diffusion and mass transport in solids, Les Editions de Physique, 1991. H. Risken, The fokker-planck equation, methods of solution and applications, 2 ed., Springer, 1996. Paul G. Shewmon, Diffusion in solids, McGraw-Hill Book Company, 1963. D. Tabor, Gases, liquids and solids and other states of matter, third ed., Cambridge University Press, 2000. Alan J. Walton, Three phases of matter, 2 ed., Oxford University Press, 1992. H. Wipf, The gorsky effect, electrotransport and the thermotransport of hydrogen in metals, Journal of Less-Common Metals 49 (1976), 291–307. Dennis G. Zill and Michael R. Cullen, Differential equations with boundary-value problems, International Thomson Publishing Inc., 1997.


Shared By:
Description: Numerical Solution to the Drift-Diffusion Equation of Point