ALGORITHMS FOR MULTI-DIMENSIONAL ANALYSIS OF SEMICONDUCTORS DERIVED FROM FIRST
PRINCIPLES
Zvonko Fazarinc
KEYWORDS: semiconductors, semiconductor structure, computer analysis, algorithms, C-language, charge carriers, transport equation, electron distribution, hole distribution, electric field, basic principles, multi-dimensional
ABSTRACT: Expressions for multi-dimensional analysis of semiconductor structiires in the discrete domain are derived from first principles. A simple structure is analyzed and the algorithms are cast in C-language.
Algoritmi za večdimenzionalno računalniško analizo polprevodnikov iz osnovnih fizikalnih principov
KLJUČNE BESEDE: polprevodniki, strukture polprevodniške, računalniška analiza, algoritmi, C jezik, nosilci nabojev, enačba transportna, porazdelitev elektronov, porazdelitev vrzeli, polje električno, principi osnovni, večdimenzionalnost
POVZETEK: Izrazi za večdimenzionalno računalniško analizo polprevodniških struktur so izvedeni iz osnovnih fizikalnih zakonov. Uporaba je ilustrirana na enostavnem primeru in algoritmi so prikazani v C-jeziku.
1. Introduction
Multi-dimensional analysis of semiconductor structures is commonly deferred to pre-canned computer programs [1] which often drape a veil of mystery over the inner workings of such design tools. The Poisson's equation is usually taken as the basis for evaluation of Fermi levels which then control the distribution of charged carriers. When the transport equation is used instead, the carrier distributions are computed from its discrete counterpart. This is prone to producing wrong answeres and makes the imposition of boundary conditions quite difficult.
We will be making use of the classical physics principles which are applicable to semiconductor structures larger than a few tenths of microns. For smaller structures quantum mechanics must be invoked and the reader should take note of this.
2. The Method of Approach
Before we address the general case we introduce the methodology with a simple example. Fig. 1 illustrates three points in space separated by Ax. With each location we associate a particle count C(x,t) at time t.
The purpose of this paper is to derive the relevant equations for multi-dimensional analysis of semiconductors in discrete form directly from first physical principles. Such approach effectively avoids the hazards of discretization of partial differential equations [2], makes the imposition of boundary conditions, intuitive and, most importantly, it provides the practicing engineer and the novice with an insight which enables them to independently access the analytical powers of computers. The limitations of a paper prevent us from developing anything resembling a complete source code. Nevertheless, we will address the crucial ideas and make them understandable so that they can be embellished with refinements when needed.
C{x-Ax,t)
x-Ax
C(x,t)
x -5>-
Cix + Ax,t)
x + Ax
Fig. 1: Illustration of three points in one-dimensional space.
Particles are assumed to be in random thermal agitation which implies that they are equally likely to move to the left or to the right. We denote the likelyhood of their motion in one and the other direction by I. This means
that /-times the number of particles in a given position will move to the left and the same number will move to the right. We make a simplifying assumption that particles which do move make one single space step Ax in one time increment At. Armed with this information we can entertain the following question: "What will the concentration be in position x at time f-t-Af given the status at time /?" The answer proceeds along the following line of reasoning. If the likely hood of motion is / then /4imes the number of particles in position x - Ax will move into position X during one time interval At. During the same time /-times the number of particles will move into position X form X + Ax. 2/-times the number of particles initially residing at xwill have moved out of this position. What we have then left at x is
C(x,/ + AO = C{x,t)+lC(x-Ax,t) +lC{x+Ax,t)-2lC{x,t)	(1)
The mathematical manipulation below is intended to show that (1) is the discrete form of the diffusion equation with diffusivity D given by
D = lAx^/At
(2)
First we subtract C(x,t)or\ both sides of equation (1) and divide by At. Then we multiply and divide the RHS of the resulting equation by Ax^and obtain
C(x,t+At)-C(x,t) ^ At
lAx^ C(x-Ax,t) + C(x + Ax,t)-2C(x,t) At	Ax^
We recognize the numerator of the LHS of the above equation as the temporal difference of C(x,t) and the numerator of the RHS as the second spatial difference of C(x,t). Upon taking the limit as Ax and At go to zero we end up with the familiar basic diffusion equation for which the relationship (2) applies.
d t	dx^
A temporal sequence of plots produced by (1) when the initial distribution is a 5-function in the center and a unit step at the left is shown in Fig.2,
An implementation of our example in C-language is shown below
// Initialization:
for (X = 0; x<200; x+ -i-) C[x] = 0; for (x=0; x<20; x-t- -t-) C[x] = 1.0; C[100]= 10.0;
//Time loop of N passes: for(t = 0;t<N;t+ +) {// Space loop: L= 1.0;
		I'O
		^IT
1 =11		
l\\\	i	w
Fig. 2: Solution of discrete diffusion equation in one dimension
for (X = 0; x< 200; x-t- +) {temp = 0.25*C[x-1] + 0.25*C[x+1] -n 0.5*C[x]; qx-l] = L; L= temp; // Plot temp vs x
The timing required by (1) has been implemented in the above algorithm by reliance on memory management of computers. The values on the RHS of an equation are always taken to be the old values and those on the LHS of equation as the new ones. But we update C[x] by the new value only after we have computed C[x+1] which requires the old value of C[x]. This is accomplished by swapping the temp variable with L every time we pass through the space loop. It is obvious that the / value has been taken to be 0.25 in the example. One may wonder if such simple algorithm truly represents the solution of the diffusion equation. A quantitative comparison of (1) with the continuous diffusion equation is given in reference [2],
Next we allow a force F, indicated in Fig.1, to act uniformly on ail patlicles. It is not difficult to conclude that a positive force - one that pushes to the right - will increase the likelyhood of particles moving to the right. At the same time it will decrease their chances of moving to the left. We can modify equation (1) for this case by introducing a skew factor /which is related to the force Fand which biases the likelyhood I in the direction of force F.
Cix,t+At) = C{xj) 4- {I +f)C{x-Ax,t) (l-f)C{x + Ax,t) - 2lC{x,t)
(3)
The effect of force on C(x,t) cancels out in this simple example as the reader may verify. Expression (3) is the transport equation in one dimension with constant diffu-
sivity and force. To prove this we subtract C(x,t) on both sides, divide by A/and separate the /and ^contributions
- 2Cix,l)]+-^[C{x -Ax,!) - +AaV)]
Now we need a relationship between the factor /and the force F. This is best done by comparing the energies involved during the move. The thermal energy kT is related to I in the same manner as the energy derived from the force field is related to f. A formal expression to that effect is
!:f = kT:
fax
It yields the following for the force factor FAX
/ = /
2kT
(4)
We substitute (4) into our last equation, multiply and divide the first term on the right by Ax^, and the second term by 2Ax and end up with the following difference equation
C{x,i+At)-C{x,t) ^ At
lAx^ C(x-Ax,t) + C(x + Ax,t)-2C(x,t) At	Ax^
. 2/Aa:^ F C{x + Ax,t) - C(x - Ax,t) At TkT	lAx
/ - 0 4/
1 \ 1 \ .
Fig. 3:
Solution of transport equation for one-dimensionai case.
(3) has been derived make the introduction of boundary conditions particularly easy.
We will now exploit the method just introduced for answering some of the more profound questions. One of them involves the diffusivity gradient. If D is a function of space, does it appear under the first or the second derivative sign? What if the diffusivity and force have different gradients along the spatial directions? How do we handle multidimensional analysis of semiconductors in general? These and some other questions will be addressed as we proceed.
The second term on the right is readily recognized as the first central difference in x of C(x,i). It becomes the first derivative with respect to x when the limit is taken. The above equation then assumes the form
dC(x,t) d'C(x,t) F dC(x,t) d t	dx^ kT dx
In the above we have substituted (2) for /Ax^/At. Fig.3 shows a plot of expression (3) as function of x for a uniform force F with time being a parameter. Initial conditions are identical to those in Fig.2.
The numeric values for the plot of Fig.3 were generated by a source code identical to that shown earlier. The only difference is that the skew factor is introduced. We have chosen its value to be / = 0.4 / = 0.1. Consequently the only modification of the algorithm is in the space loop which now reads: temp = 0.35*C[x-1] + .15*C[x-(-1] + .5*C[x]; Everything else remains the same. The reader, familiar with problems arising in computer solutions of transport equation may find (3) to be of considerable interest. Its simplicity and the physical basis from which
3. Transport of Particles in Two Dimensions
The motion of charge carriers in semiconductors is governed by thermal energy, by electrical forces* and by properties of material through which they are moving. Their number depends on influx and outflow and on generation and recombination of oppositely charged pairs. We will derive the relevant equations taking into consideration the spatial variability of these effects and will allow, in addition, the temporal variation of electric fields and of carrier concentrations. In order to shield the derivation from excessive notational complexity we will limit it to the two-dimensional case. The extension to three dimensions will become self evident as we proceed.
We start with the illustration in Fig.4 which shows two points in x-space separated by Ax and a third point displaced in y-space by Ay. We assign the likelyhood of thermal motion in x-direction by Ix and that in y-direction
* Gravitational forces are negligible in comparison

• h. f.. • X-Ajco'	Jf^
Equations (6) and (7) contain the particle flux information which we intend to address later. At this time we remain focussed on our goal to obtain the complete transport equation of the form (3) in two dimensions and for
ly- fy. Fy.,Ny.

Fig. 4: Illustration of three points in two-dimensional space.
by ly. We also allow a force Fx to act in the x-direction and a force Fy in the y-direction. These forces are responsible for the respective skew factors fx and fy through the relationship (4).
The motion likelyhood I and its skew f are assigned to the space between discrete spatial locations rather than to the locations themselves. The reason for this can be understood physically if we considerthe case where the point (x-Ax,y) belongs to the conducting region and the point (x,y) to an insulating layer characterized by / = 0. The particles can move into the insulator but they could never escape if the property of space /= 0 were assigned to (x,y) itself. On the other hand, if / = 0 is assigned to the space between the two points, no particle exchange can take place across the boundary. Force factor f, being a bias to I must share the same space. The rule is implemented in Fig.4 by making the subscript "x" to mean "half Ax to the left of x" and "y" to mean "half Ay below y. Using this notation we can state that during one time interval Af there will be [lx(x-,y) + fx(x.,y)]C(x-Ax,y,t) particles flowing from left to right and [lx(x. ,y)-fx(x-,y)]C(x,y,t) particles flowing from right to left. The net numberof particles /Vxfx.,yjf lowing from (x-Ax,y)\o (x,y) is the difference of these two terms
x^+Ay
^x{x.,y) = Ux.,y)+fAx.,y)]Cix-Ax,yj) -UAx.,y)-fAx.,y)]C(x,yj)
(6)
Similarly we get for the vertical or y-direction the following expression for the net number of particles Ny(x,y-) flowing from (x,y- Ay) towards (x,y)
Ny(x,y.) = [ly(x,y.)+fy{x,y.)]C(x,y-Ay,t) -[lyix,y.)-fy(x,y.)]C(x,y,t)
(7)
^y + > +
fy*
•	fx. •
Fig. 5: The other points in two-dimensional space
variable diffusivities and forces. This requires that we know all the influxes and outflows from (x,y). Fig.4 guided us to obtain two of them, expressed in (6) and (7). With the aid of Fig.5 we obtain the other two which flow to the right and upward of (x,y) respectively.
We denote the net number of particles flowing out of (x,y) towards (x + Ax,y) by A/xfx+.yj. Their count is
NAx + ,y) = [I, ix + ,y) (x + .y )]C {x,y,t)
- iUx + ,>')■ fx(x ^ ,y)]Cixi-Ax,y,t)	^^^
The net number of particles Ny(x,y+), flowing from (x,y) upward towards (x,y+ Ay) is
Ny{x,y + ) = [ly(x,y^)+fy(x,y^)]Cix,y,t)
-[ly(x,y,)-fyix,y,)]Cix,y+Ay,t)
We can write an expression for the particle count at time f + Af at the point (x,y) by adding all inflowing particles to C(x,y,t) and subtracting all outflowing particles from it. We will allow for the possibility that in addition to particle exchange some of them are being generated at the rate G(x,y) while some are being lost at the rate R(x,y). This produces a net particle increase of [G(x,y) -R(x,y)] At during the time interval At. The total particle count at (x,y) at t+ At. The total particle count at (x,y) at t+ At'is then
Cix,y,t+At) = C{x,yj)
+NAx.,y) + Nyix,y.)-NAx^,y)-Myix,y^)
-^[G{x,y)-R{x,y)]At
Equation (10) in association with (6) through (9) is the discrete form of the two dimensional transport equation and we will make use of it shortly. But first we will perform a limiting process on it to obtain the equivalent differential equation. Rather than blindly substituting expressions (6) through (9) into (10) by noting that the difference Nx(x+,y,t) - Nx(x-,y,t) refers to the midpoint between (x.,y) and (x^,y). But this is excatly (x,y) by our definition of x-and Consequently we can write the difference Nx(x^,y,t) - Nx(x-,y,t) as ^xNx(x,y,t)^N\^ere Ax signifies that the difference is with respect to x. When we do the same with the y-difference we get for (10)

(11)
The value for Nx(x,y,t) can be formally derived from (6) or (8) by respectively incrementing or decrementing all x-arguments by x/2. The result is
C{x^,t +A0 = C{x,y,t) + - Ay[//A:j')2C (rj-.Ol + [G (^J-)
By means of (4) we convert /-factors to /and then employ (2) to convert /-factors to diffusivity D. We also subtract C(x,y,t) on both sides of equation and divide them by At. A straightforward algebraic manipulation leads to
C(x,y,t+A().C(x.v.t) ^ ^ Ar	Ax
A.
Aa:

Dxix^) A,


Ay
kT
^vC-^iV.O
Ay
Dyix^)
Ax
AyC(X^,t)
Ay


(12)
In the limit when all increments go to zero the above becomes a differential equation
In the above we have used the established notation for half-Ax values as, for example, C(x-Ax/2,y,f) = C(x.,y,f).
Similarly we get by either incrementing (7) or decrementing (9) by Ay/2 the expression for Ny(x,y)
Ny(x,y,t) = [ly{x,y)+fy(x,y)]C(x,y.,t) -[ly{x,y)-fyix,y)]C(x,y^,t)
(13)
Substitution of (12) and (13) into (11) yields
-Ax fAx^)lCix.^,t) + Cix,o',t)]
-Ay ly(xo>)[C(x,y.,t)

■R(x^)]At
This time we recognize the term C(x-,y,t) - C(x,,y,t) to be the negative difference centered on (x,/). Consequently we can denote it by - AxC(x,y,t). Similarly for the y-direc-tion. Furthermore the sum of C(x.,y,t) + C(x^,yt) = 2C(x,y,f) and becomes exact when the increments go to zero which we are just about to do. We have now
9C(MA1 = Aß
dl
ate
{x,y)

ax
kT
IH

{j'^x)j(<cxyj 'Wxfoe + G{x,y)-R{x,y)
+
While this equation may be the correct starting point for a computerized solution of the transport equation when diffusivities, forces and concentrations vary in both dimensions it certainly does not guarantee the correctness of the solution. There are countless possible ways to discretize a differential equation but very few of them yield correct answers. Therefore it is strongly advisable to start from equation (10) which is derived from fundamental physical principles and circumvents the discretization problem altogether.
4. Flux of Particles in Two Dimensions.
Equations (12) and (13) provide us with the number of particles moving in the x and y direction, respectively, across the point (x,y). The corresponding flux is that number multiplied by the particle velocity. Our initial supposition was that particles move one space interval in one time interval At. Consequently their velocity is Ax/ At in the x-direction and Ay/ At in the y-direction. The two respective fluxes are then


At	Ax

The same kind of reasoning which led to the last differential equation in the previous Section guided us from here to the particle flux equations in differential form

= -Dyix^)

If we multiply the flux by electric charge q we get the electric current density. The forces in this case are qEx and qEy, respectivelyforpositive charges, i.e., holes and -qEx and -qEy for electrons where Ex and Ey are the electric fields. The concentrations C(x,y,t) are n(x,y,t)ior electrons and p(x,y,t)io{ holes. The corresponding four current density equations are
inxixo',0 = qD^{x,y)
ipx{x,y,t) = -qDp,{x,y) }py(x,y,t) = -qDpy(x,y)
dx kT
ctx kT
The reader may note that the so called Einstein relationship between carrier diffusivity D and their mobility |a is contained in the above equations
= D "
kT
Before we leave this section let us examine equations (12) and (13) for physical consistency. Assume that the force factor f which depends on the electric field becomes larger than / in (12) or (13). This implies that the particle flow represented by the second term of either equation reverses its direction. This could quickly deplete the particles at the adjacent point and would even make their concentration go negative. Such physical impossibilities usually present themselves as numeric
instabilities. We can therefore derive a stability condition from requiring that f never exceedes the value of I. Using (4) this condition can be writen as FAx/kT< 1. For the electric case this translates into EAx<kT/q. Because E Ax represents the potential difference over one grid point our stability condition requires that this potential difference never exceeds 25 mV at room temperature.
5. Forces on Particles in Two Dimensions.
The last unresolved quantity that appears in equations (6) through (9) is the force factor f. According to (4) this is directly related to the electric field for the two respective directions. The electric field, on the other hand, emanates from electric charges as suggested by the Coulomb's law. Before we apply it to our case a discussion is in order.
In real, three-dimensional space the electric field surrounding a charge possesses spherical symmetry and is as such decaying as the square of the distance from the charge. If we were dealing with a three-dimensional case we would compute the field at a given position (x,y) by adding vectorially the contributions from all surrounding charges using an inverse square law. Then we would resolve that field into its Cartesian or other components as needed by the analysis. Instead of spherical symmetry we have a cylindrical symmetry for our two-dimensional case and consequently we can allow the electric field to decay only linearly with distance away from the source charge. This may appear unconventional but is no more so than is a two-dimensional space. In a one-dimensional space the field does not decay at all and we therefore compute the field as the sum of net charges without prorating their effects for distance. One must be quite cautious when applying natural laws to unnaturals spaces. With this in mind we continue now with the evaluation of force factors to be applied to equations (6) through (9).

/

• / •
{x-Ax.y) (x.,y) (x.y)

(.x.y-C^y)
Fig. 6: Aid to calculation of the electric field components
We remember that the force factors have been defined halfway between discrete points. With the aid of Fig.6 we compute the x-component of the electric field at the point (x,y) as produced by the charge Q(u,v) a distance
rx away. A linear decay of the field with distance is ennployed as appropriate for our two-dimensional case.
6. An Example of a Two-Dimensional Semiconductor Structure
E,{x.,y) = E(x.,y)
Q(u,v) "-a:.
/V	n
u -X.
QM -
2jr	+
Similarly we get for the y-component at (x.y.) according to Fig.6

I'KTy	Ty
v-y.
QM -
27r {v-y.f^{u-xf
The total field is the summation of components due to all charges Q(u,v)
EÄx.,y) = EEEE
2n	2 ^ (16)

and
Ey{x,y.) = EEEE

We will now illustrate how (10) can be applied to impurity redistribution and to subsequent charge migration and how (16) and (17) can be used to compute the resulting forces opposing such migration. This is done in the hopes that the reader will be encouraged to experiment with multidimensional electric transport problems on his own.

I C(A-,y, 5mm)
We start with the doping profile CCx.yj shown on the LHS of Fig.7. It represents a PIN junction with a very narrow intrinsic region, with a donor density of Wd=4*10 and an acceptor density of Na=2* 10^®. The profile shows the difference C(x,y) = Noix.y) - na(x,y).
The widths of the two regions are approximately 2.5|i each. Next we anneal the sample by subjecting it to 1100°C for about 5 minutes. The diffusivity was chosen to be 5* 1at this temperature which is representative of commonly used dopants. The redistributed impurities are shown on the RHS of Fig.7.
This result was achieved by using (6) through (10) with the following parameters:
Ix = I y = . 125 independently of x and y, and fx = fy = 0 and G - R = 0 everywhere.
Equation (10) simplifies then into
C{x,y,t +At) = 0.5C(x,>',r) + 0.125[C(j:-Ax,.')
^C{x+Ax,l) + C{x^- Ay,t) + C(x^+ Ay.O]
(18)
Fig. 7: Impurity prof He before and after annealing.
Denote the number of Ax steps by X and the number of Ay steps by Y and assign a two-dimensional array C(x,y) to the impurity concentration. We can then write the algorithm which represents the above equation in C-lan-guage as
for (y = 0; y<Y; y-n-) L[y] = 0.0; for (x = 1; x < X; x++) (for (y = 0; y < Y; y+-i-) {temp = .5'C[x,y] + . 125*(C[x - 1 ,y] + C[x + 1 ,y] -I- C[x,y - 1] + C[x,y + 1]);
C[x- 1,y] = L[y];L[yJ = temp; }
} for (y=0; y < Y; y++)
{C[X -1,y]= L[y];L[y]= 0.00;}
A few explanations are in order. First, the choice of I seems fairly open since it represents the diffusivity in conjuction with space and time steps according to (2). But if you consider the case where all the matter is initially concentrated in position (x,y) then we can have no more than 1/8 of the central matter moving out in each direction while preserving an equal amount of matter in the center. If we allow higher depletion, instabilities may occur. The rule of thumb is then to choose /max=0.125 for highest diffusivity and proportionately less for all other values of D. In a three-dimensional case the maximum value of /should be chosen to be 1/12 and for the one- dimensional analysis /ma;<=.25.
For our example we have chosen 25 segments along the entire width and 25 along the entire lenght of the
sample corresponding to a space step of 0.2|a . According to (2) the time step At for our chosen diffusivity amounts to 100 seconds, meaning that every time we traverse the algorithm the time advances by this amount. Consequently, the annealed profile in Fig.7 is achieved by executing the algorithm only three times (five minutes).
In the C-algorithm above we have taken care of timing implied by (18) which demands that all RHS variables be taken at time t. As the execution moves from lower to higher values of x and y we must postpone the update of concentrations until we have moved to the next x-col-umn. Therefore we save one whole x-column in L[y] and load it into C[x -1 ,y] column after the old values are no longer needed. The first "for-loop" initializes the L-array to whatever boundary value we wish to assign to the x = 0 column. As the x-loop starts with the index 1, the column x = 0 receives the initial values contained in L[y]. The very last "for-loop" updates the last x-column at X-1 and resets the L-array to the initial value which will subsequently restore the boundary value at the x = 0 column again. Whatever we load into the X column which lies outside the reach of the x-ioop becomes the upper boundary value for x. This column never gets updated so it supplies the boundary conditions which it contains to the X-1 column. In a similar fashion we can impose the boundary values to the y-column or to any point in the x-y domain, for that matter. This we do by loading such conditions into the appropriate location everytime the loops are traversed and the initial value there has been changed by the algorithm.
This lenghty discussion of the boundary value problem is intended to encourage those who have experienced difficulties with similar problems in the continuous domain. The discrete nature of computerized analysis makes the imposition of boundary conditions almost intuitive.
Next we must compute the concentrations of holes and electrons produced by the annealed profile. This is quite straightforward when we invoke the neutrality condition Nd-Na +p-n = 0 and the mass- action law np = n^[3,4],
The following algorithm has been used to obtain Fig.8

p{'t.y)
for (x = 0; X < X; x++) {for(y= 0;y<Y;y++) {if (C[x,y] > 0.0)
{n[x,y] = C[x,y]/2,0 + sqrt(C[x,y]*C[x,y]/4.0 2.0e20);
p[x,y] = 2.0e20/n[x,y]; }
if(C[x,y]< 0.0)
{p[x,y] = - C[x,y]/2.0 + sqrt(C[x,y]*C[x,y]/4.0 + 2.0e20);
n[x,y]= 2.0e20/p(x,y); }
if (C[x.y]= = 0.0) (n[x,y] = 0.0; p[x,y] = 0.0;} }
}
where 2.0e20 stands for squared intrinsic concentration of Silicon.
Gradients of holes and electrons seen in Fig.8 give rise to their migration into the adjacent regions of low density which destroys the initially imposed neutrality. This gives rise to electric fields which tend to oppose the migration and which we will compute later. First we use (10) again to compute the movement of charges but this time we want to distinguish the migration rates of holes and
Q{x,y) = C{x,y)+p{x,y)-n{x.y)
Fig. 9: Net charge distribution after 10 picoseconds
electrons. We assign a diffusivity Dn = 28 to the latter and Dp = 9 to the former. These numbers are in fair agreement with the respective impurity concentrations. We assign the value ln= .125 to the highest diffusivity. The corresponding value of Ip must then be .125x9/28 = .04. With these values and with the force factors initially set to zero we have obtained the net charge density distribution Q(x,y) shown in Fig.9 after approximately 10 picoseconds.
It is apparent from Fig.9 that the departure of electrons along the four boundaries of the A/-region produces a net positive charge. Similarly the departure of holes from the P-region produces a net negative charge all around the boundary. A charge reversal is observed at the junction. Implemention in C-language of equations (8) through (10) which have produced Fig.9 is shown below.
Fig. 8: Electron and hole distribution.
for (i = 0; i < 5; i++) {for (x= 0; X < X; x++) {for(y= 0;y<Y;y-f-+) {Fxy = .125*(fx[x,y] - fx[x+1 ,y] + fy[x,y]
-fy[x,y+i]);
tempN = (.5-Fxy)*n[x,y] -f- .125*(1.0-fx[x,y])*nlx-1,y] + .125*(1.0 + fx[x-f-1,y])*n[x + 1,y] + .125*(1.0-fy[x,y])*n[x,y-1] + .125*(1.0 + fy[x,y+1])*n[x,y+1]; n[x-1,y]= L[y];L[y]= tempN; tempP= (.84 -t- Fxy)*p[x,y] -t- .04*(1.0 + fx[x,y] *p[x-1,y] -t- .04*(1,0 - fx[x+1,y]*p[x+1,y] + .04*(1.0 + fy[x,y])*p[x,y-1] + .04*(1.0 - fy[x,y+1])*p[x,y+1]; p[x-1,y]= L[y];L[y]= tempP;
}
} }for(x= 0;x< X-,x++) for (y= 0; y < Y; y++) Q[x,y] = C[x,y] + p[x,y] - nlx,y];
In the above we have made the following substitutions fx[x,y] = E,{x.,y)Ax^ bc[x+ l.y] = E,{x ,

kT
kT

delX= u-x; if(delX!= 0) {delY= v-y + .5;
EY = EY + Q[u,v] *delY/ (*delY*delY -)- de!X*delX), }
}
fx[x,y]= EXM.Oe-10;
fy[x,y]= EY*1.0e-10; }
The factor 10"^° in the above takes account of the permittivity of Silicon, of the integration intervals of the factor q/kTanö of the conversion between surface and volume density of charges involved. An interpretation of the major features of Fig. 10 tells us that the force on electrons and holes is directed inward both along the x-axis. The y field shows more intensity in the Af-region than in the P which has its origin in the higher donor concentration and electron diffusivity. One major but expected feature of the electric field Ex is its high intensity along the junction region. This field eventually blocks the further migration of electrons and holes into the adjacent regions. An actual simulation of a semiconductor structure would alternate between the last two algorithms, computing the net charges from the first one and finding the resulting force factors from the latter one. These would then be substituted back into the first one. An equilibrium would eventually be established at which point te migration would cease completely.
E.U^y)
Fig. 10: Components of the electric field.
We have allowed five passes through the algorithm. Because the time interval is Af = / CAx//D each pass is somewhat shorter than 2 picoseconds. The net charge density Q(x,y) is obtained from the resulting distributions of p(x,y) and n(x,y) as the last step in the algorithm. When we substitute this into (16) and (17) we obtain the two components of the electric field. These are shown in Fig. 10 and a C-version of the algorithm which generated the numeric values is reproduced below the figure.
for (y = 0; y < Y; y++) for (x = 0; X < X; x-n-) {EX= EY= 0.0; for (u = 0; u < X; u++) for (v = 0; V < Y; v+-t-) {delY= v-y; if(delY!= 0) jdelX = u - X + .5;
EX = EX -H Q [u,v] *delX/(delY*delY + delX*delX); }
Before we conclude let us recognize that the electronic fields could be computed in an alternate way. It is known [5], [6], [7] that the expression we have derived for the electric field from Coulomb's law is in fact the formal solution of the Poisson's differential equation
ax dy . e
Where Q is the charge density. Unfortunately this one equation with two unknowns and as such not very useful. But it is also known from the field vector theory [8] that field components Ex and Ey are the gradients of a function V(x,y) which satisfies the Poisson equation in the following way
dx^	dy^	e
When this equation is solved for the scalar quantity V(x,y) the field components can be obtained as the respective partial derivatives

and
Ey(X^)
.SEML
dy
It turns out that the effort in solving the Poisson equation for the potential is comparable to what we have done and our approach is more in keeping with the promise of staying close to first principles. Nevertheless it ought
to be pointed out that the Poisson equation in terms of potential V(x,y) can be solved by the same algorithm we have developed for the particle transport by diffusion. The diffusion equation in two dimensions with constant diffusivity has the differential form
dV{x,y,t) d'V(x,y,t) , d'V(x,y,t) dt	dx^	dy^
When add a source function ^^ to the RHS of the
above and solve the equation for steady state condition, i.e. dV(x,y,t)/Qt{he equation goes over into the form
taking advantage of present-day computer performance.
Ackknowledgement
The highly efficient and versatile algorithm which was used to plot the two-dimensional distributions was developed and given to me by Ing. Dragan Fidler of Zemun.
References
d'V{x,y) , d'V(x,y) Q{x,y)
e

Consequently we can solve the Poisson equation with the aidof equation (15) as implemented inthe respective C-algorithm when we add the source function Q(x,y)k to the RHS of the temp expression. When there are no more changes between two successive evaluations of temp the resulting distribution C[x,y] is the solution of the Poisson equation. The results obtained for our test case are identical with one or the other method and the time consumed is about the same.
7. Conclusions
The transport equation and the equivalent of the Poisson equation were derived from first principles without invocation of unsubstantiated abstraction. The two-dimen-sional case has been tested on a simple semiconductor structure and the corresponding algorithms have been presented in the C-language. The intent was to show that the availability of computers warrants a fresh look at the traditionally accepted mathematical models and that more intuitive approaches are made possible by
[1]	SUPREME by Stanord, PISCES and DEPICT by Tehnology Modeling Associates, SIMPEL by Berkeley, MASTER by Silvacor.
[2]	Z,Fazarinc, Discretization of Partial Differential Equations for Computer Evaluation, Computer Applications in Engineering Education. Vol.1, No 1, pp.73-85, John Wiley, 1992-93.
[3]	0.Kittel and H.Kroemer, Thermal Physics 2-nd Ed., Freeman & Co., San Francisco 1980, p270 Derived from Thermodynamics, p268 from reaction rates, p361 -3 from density of states.
[4]	A.S Grove, Physics and Technology of Semiconductor Devices, John Wiley 1967, p 100-1.
[5]	James 0. Maxwell, A Treatise on Electricity & Magnetism, Clarendon Press, 1891, Chapter II.
[6j Oliver Heaviside, Electrician, Phil. Mag. XXVII, 1889, p.324
[7]	Simeon D.Poisson, Bull.de la Soc. Philomathique iii, 1813, p.388
[8]	I.Sokolnikoff and R.Redheffer, Mathematics of Physics and Modem Engineering, McGraw Hill, 1958, p386-7
prof. Zvonko Fazarinc, Ph. D Stanford University, ret. IHewlett-Packard Laboratories, ret.
880 La Mesa Drive Menio Parl<, CA 94025
Prispelo: 27.08.93
Sprejeto:17.09.93