Diffusion 

Ping Du
Manufacturing Engineering
Boston University
duping812@gmail.com 

June 27, 2007 

Diffusion is the movement of particles from an area of high concentration to an area of lower concentration until equilibrium is reached. It can be modeled quantitatively using the diffusion equation, which has different forms in different physical situation. Such as molecular diffusion is governed by Fick's law, thermal diffusion (aka thermal conduction) is governed by Fourier's Law, and electron diffusion in an electrical field essentially turns to Ohm's law[1]. In this paper, doping process in semiconductor industry is chosen as an example to illustrate some  basic ideas of the physics and understand Maple's capability to solve partial differential equations (PDE). 

Note: This paper is also presented as the final project for the course MN 526, Simulation of Physical Processes. 

Initialization 

>	restart: with(plots):

 

>	unprotect(D);

 

>	alias(c[0]=c0, c[1]=c1, c[2]=c2);

 

(1.1)

 

1. Diffusion equation 

1) Fick's first law 

It states that the flux of the diffusing material in any part of the system is proportional to the local concentration gradient:  

 

where is the flux, D is the diffusion coefficient, is the concentration. 

 

2) Fick's second law 

It states that a change in concentration in any part of the system is due to inflow and outflow of material into and out of that part of the system. Effectively, no material is created or destroyed:  

 

if D is a constant, then the equation can be simplified as linear:  

= 

Consider the simplest 1-D condition:  

>	PDE:=diff(C(x,t),t)=D*diff(C(x,t),x,x);

 

(2.1)

 

2. Analytical solutions 

In semiconductor industry, doping  is a typical diffusion process, either by means of ion implantation or vapor phase deposition[2]. By simplifying the boundary conditions for these two procedures, the particular analytical solutions of diffusion equation can be derived. There are many mathematical ways to derive them, such as Fourier transform, which could also be done in Maple[3]. 

 

1) Gaussian solution 

First implant a very narrow peak of dopant, which approximates a delta function of concentration. Taking the origin to be at the delta function, the boundary conditions are: 

 

 

and 

 

where Q is the total dose contained in the peak. Q is fixed and remains at all times. The solution satisfying these conditions is:  

>	f1:=Q/(2*sqrt(Pi*D*t))*exp(-x^2/(4*D*t));

 

(3.1)

 

>	c1:=unapply(f1,x,t);

 

(3.2)

 

This solution can be checked by taking derivatives in both sides. 

>	Diff(c1(x,t),t)=diff(c1(x,t),t);

 

(3.3)

 

>	combine(simplify(rhs(%), radical));

 

(3.4)

 

>	D*Diff(c1(x,t),x,x)=D*diff(c1(x,t),x,x);

 

(3.5)

 

>	combine(simplify(rhs(%), radical));

 

(3.6)

 

From above results, the two sides are identical. 

 

2) Error-function solution 

The other solution is the situation of diffusing from infinite source of dopant, which might correspond to vapor phase deposition. In this case, the boundary conditions are: 

 

 

This situation can be considered as a sum of a series of slices, each of width Δx and concentration , and satisfying the previous Gaussian solution. Using linear superposition and taking the limit of Δx, the sum becomes an integral, which is also called the error-function: 

>	f2:=c0/2*erfc(x/(2*sqrt(D*t)));

 

(3.7)

 

>	c2:=unapply(f2,x,t);

 

(3.8)

 

>	Diff(c2(x,t),t)=diff(c2(x,t),t);

 

(3.9)

 

>	combine(simplify(rhs(%), radical));

 

(3.10)

 

>	D*Diff(c2(x,t),x,x)=D*diff(c2(x,t),x,x);

 

(3.11)

 

>	combine(simplify(rhs(%), radical));

 

(3.12)

 

3. Analysis 

1) Gaussian solution 

There are several immediate consequences of this solution. The peak concentration locates at x=0, and decreases as 1/ 

>

c1(0,t);

 

(4.1)

 

When the distance from the origin x=2the concentration has fallen by 1/e. 

>	c1(2*sqrt(D*t),t)/c1(0,t);

 

(4.2)

 

Actually the above two quantities are very similar to the mean and standard deviation in normal distribution. 

 

Now substitute some values for the coefficients to check the resulting concentration profile. 

>	c0:=1; Q:=1; D:=1;

 

 

 

	(4.3)

 

>	plot([c1(x,0.1), c1(x,1), c1(x,10), c1(x,100), c1(x,1000)], x=-10..10, legend=[0.1, 1, 10, 100, 1000]);

 

 

2) Error-function solution 

>	plot([c2(x,0.1), c2(x,1), c2(x,10), c2(x,100), c2(x,1000)], x=-10..10, legend=[0.1, 1, 10, 100, 1000]);

 

 

4. Numerical solution 

1) First try solve the PDE using the build-in function pdsolve. 

>

PDE;

 

(5.1)

 

>	pdsolve(PDE);

 

(5.2)

 

2) It shows that Maple could not find the exact analytical solution, then try solve it numerically.  

Consider a simple situation, take the boundary condition to be left end fixed at constant concentration and right end insulated. 

>	IBC:={C(x,0)=0, C(0,t)=c0, D[1](C)(10,t)=0};

 

(5.3)

 

>	pds:=pdsolve(PDE,IBC,numeric);

 

(5.4)

 

>	L1:=[0.01, 0.1, 1, 5, 10]; L2:=[red, green, yellow, blue, magenta, black];

 

 

	(5.5)

 

>	for i from 1 to 5 do  pn[i] := pds:-plot(t=L1[i], color=L2[i]): end do:

 

>	display({seq(pn[i], i=1..5)}, title=`Numerical solution at t=0.01, 0.1, 1, 5, 10`);

 

 

3) Compare with analytical solution.  

The above boundary condition satisfies the error-function solution near a surface. Due to symmetry, the midpoint in the original plot now becomes the left endpoint in the new plot. Thus doubling the original left endpoint could lead to the same solution. 

>	f2_s:=2*f2;

 

(5.6)

 

>	c2_s:=unapply(f2_s,x,t);

 

(5.7)

 

>	for i from 1 to 5 do  pa[i] := plot(c2_s(x,L1[i]), x=0..10, color=L2[i]): end do:

 

>	display({seq(pa[i], i=1..5)}, title=`Analytical solution at t=0.01, 0.1, 1, 5, 10`);

 

 

4) Now check the error between the numerical solution with the analytical solution. 

>	for i from 1 to 5 do  perr[i] := pds:-plot(abs(C(x,t)-c2_s(x,t)), t=L1[i], color=L2[i]): end do:

 

>	display({seq(perr[i], i=1..5)}, title=`Error check at t=0.01, 0.1, 1, 5, 10`);

 

 

From the above plot one can see that the numerical solution fits very well when t>1, but when t<0 the numerical solution has a large error at the very beginning. This is due to the timestep option used in the pdsolve function. Since the previous solutions at each time use the same value for timestep, the magnitude of error has more possibility to enlarge when evaluating at smaller time. 

 

The default value of timestep is same as spacestep, which is 1/20th of the space range. In this problem, the space range is 10, thus the timestep is only 10/20=1/2. That's why the solution at t=0.01 and 0.1 more error than others.  

 

5) Now try use smaller the timestep. 

>	L3:=[1/2, 1/10, 1/50, 1/100]: for i from 1 to 4 do  pds_ts[i]:=pdsolve(PDE,IBC,numeric,timestep=L3[i]):  for j from 1 to 2 do    pn_ts[i,j] := pds_ts[i]:-plot(t=L1[j], color=L2[i+j]):  end do: end do:

 

>	display({seq(pn_ts[i,1], i=1..4)}, title=`Numerical solution at t=0.01 with timestep=1/2, 1/10, 1/50, 1/100`);

 

 

>	display({seq(pn_ts[i,2], i=1..4)}, title=`Numerical solution at t=0.1 with timestep=1/2, 1/10, 1/50, 1/100`);

 

 

The above two plots show that the numerical solution at t=0.01 with timestep=1/100 approximates better to the analytical solution for satisfying this condition C(0, t)=1, while that of t=0.1 approximates earlier with timestep=1/10. It is clear that with smaller timestep, the solution is more precise. 

5. Conclusion 

In the actual situation, the diffusion coefficient D may not be an ideal constant but vary with time or position. Thus the diffusion equation turns into nonlinear form, and we can not find analytical solution easily but only use numerical methods. However, the analytical solution discussed in this paper can still serve as a guideline for further more complicated problems. 

6. References 

[1]. Diffusion, http://en.wikipedia.org/wiki/Diffusion 

[2]. James D. Plummer, Michael D. Deal and Peter B. Griffin, Silicon VLSI Technology: Fundamentals, Practice and Modeling, Prentice Hall, 2000. 

[3]. Richard H. Enns, George C. McGuire, Nonlinear Physics with MAPLE for Scientists and Engineers, 2nd edition, Birkh?user, 2000. 

Legal Notice: The copyright for this application is owned by the author(s). Neither Maplesoft nor the author are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the author for permission if you wish to use this application in for-profit activities.