Semiconductor Physics P-N junction

Thread Starter

DJ Veni

Joined Jan 9, 2007
Hi all,

I need to find a way to go about how to simulate the behavior of carrier concentration in the depletion region of a p-n junction silicon diode given the doping concentration under "non-depletion region approximation conditions" using MATLAB and thus plot the energy band diagrams and then apply it to a BJT.

Firstly Im confused as to which road to take , uniformly doped P-N where the Nd and Na dopants remain constant (FIG 1) or whether to vary them linearly while approching the junction at which Na (acceptor concentration/cm3)= Nd(donor concentration/cm3) (FIG2)---> i.e P type formed on N type substrate.



Joined Nov 17, 2003
If you could post up a diagram of what you are trying to say it would be easier to comment. What I would say is, if to create an "abrupt" P-N junction there will be some level of diffusion across the junction. If you consider one of the sides of the junction being initially doped with a particular particle concentration, then the effect of diffusion will result in a concentration gradient -dn/dt in the direction of the particle movement (i.e. the diffusion direction across the P-N junction) which closely approximates a straight line.

As for you energy bands, from what I recall the changes in the energy band characteristics will only occur over the depletion region width. I do have some work on this so should try and dig it out.


Thread Starter

DJ Veni

Joined Jan 9, 2007
I have uploaded the majority and minority carrier concentration of p and n and their concentration gradient. I have to "not use the depletion region approximation" thus im assuming that the Poissons equation remains as

p=q( Nd-Na+p-n) from which i can obtain a plot of the electric field intensity and then the electrostatic potential plot and with that the energy band diagram.
The question is how can i do it in MATLAB??
Then eventually using the drift and diffusion current quations i can make the changes to what i got.

At equilibrium
drift, diffusion current =0

Mp= mobility of holes

similarly as in the case of electrons.
I get the theory bit, but i dont know how i can implement it using matlab.

Thanks for your prompt reply dave.



Joined Nov 17, 2003
Sorry for being a touch late getting back to you.

Ok, so do you have a Mathematical representation of the Electric Field strength as a function of the distance?

Firstly, you need to define your distance vector (x) at a specified interval:

Rich (BB code):
x = [0:1e-6:30e-3];
Will generate a distance vector from 0 to 30mm (this maybe too big so you can alter your distance) at an interval of 1um (this can also be altered), giving you 30001 sampling points

Then you need to express your Electric Field equations as a function of the distance vector. As an example, if e = 2*(Sin(fx)), where f is an arbitrary variable equal to 1000, then:

Rich (BB code):
e = 2*(sin(1000*x));
Will generate the solution to e over the sample range, hence e will have 30001 sampling points.

To view this expression with the distance vector on the x-axis and the function e on the y-axis, type:

Rich (BB code):
If you provide a copy of your mathematical theory I can give you more pointers about the Matlab programming.

Since this moving more onto Matlab programming, I will move this thread to the Programmer's Corner.


Thread Starter

DJ Veni

Joined Jan 9, 2007
n= number of electrons/cm^3
p=number of holes/cm^3

under equilibrium conditions

n=p=ni= 10^(10) /cm^3 for Si-Silicon

leading to n=ni*exp((Ef-Ei)/(k*T)

where Ei= Intrinsic energy
Ef= Fermi level
k- boltzman constant
T= absolute temperature

under equilibrium:
n*p= (ni)^2

we have the charge neutrality relationship as:

Diffusion current is given as
J p/diff =-q*Dp*delta(p)
J n/diff = q* Dn*delta(n)

The total net carrier currents in a semiconductor arise as the combined result of the drift and diffusion currents. Summing the respective n and p segments we get

Jp= J p/diff + J p/drift
Jn= J n/diff + J n/drift


remember: Ef remains constant even if doping concentration varies

Under equilibrium conditions the total current remains Zero!! i.e the drift and diffusion currents will vanish only if E=0 and delta(n) and delta(p)=0------------- see equations for drift and diffusion currents

now each and every type of carrier action drift , diffusion gives rise to a change in the carrier concentrations with "Time". thus

dn/dt=dn/dt (drift) + dn/dt (diff)
dp/dt=dp/dt (drift) + dp/dt (diff)

"with this above conditions i am hoping with MATLAB i can see my plots of the carrier concentration to vary with some given time steps BUT:confused: i dont know how to use differential equations in MATLAB"

" I am assuming the Minority carrier diffusion equations will also play somepart in this"


My first step will have to get the n and p concentrations both minority and majority carrier concentrations with out the "influence of time" BEFORE i use all that previous theory to bring in the time factor and make my changes.

for my first step I need to use my Poissions equation which

Poissons equation :

Charge Density :

charge density = q(p-n+Nd-Na)

Simplified to one dimension by
the depletion approximation:

which is what i have to change from " one dimension by depletion approximation" TO " one dimesion WITHOUT the depletion approximation"

thus n and p remains even in my electric field expression!!

<<< so i need to find and in MATLAB the carrier concentration without the time factor directly and then with the poissons equation to get
plots for:

Electric Field
and as we know the inverse of the voltage curve we can get the energy band curve.

and using MATLAB i have to simulate the above theory and obtain relavent plots for my carrier concentration.


Joined Nov 17, 2003
Apologies for the lateness of my reply, I've been a little snowed in.

First thing to do is create a M-file to write your script. Open Matlab, go to File>New>M-File. Save it with a suitable name.

Declare all variables, i.e. Na, Nd, k, T, q. Your syntax should be of the form:

Rich (BB code):
% Lines preceded with one of these => is a comment
% Set ni equal to 1x10^10
ni = 1e+10;
Remember to put the ; at the end of the line to suppress Matlab returning the result to the Matlab command line.

Calculate the simple equations using your above variables:

Rich (BB code):
p = ni*(exp((Ei - Ef)/(k*T)));
As for differential equations, you state that "you to get the n and p concentrations both minority and majority carrier concentrations with out the "influence of time" - I cannot help you with this because I don't have a grasp of what you are trying to achieve here. The way I see it, you need to express n in terms of a time variable t. Define t over some length:

Rich (BB code):
% Define t from 0 to 60 second at intervals of 0.5 seconds
t = [0:0.5:60];
You will then need to ascertain what n is at time 0, 0.5, 1, 1.5 etc. Though from your description I cannot ascertain what that is. ultimately you will have a vector, n, which is 1x121 points long - your time representation of t as a vector (Matlab's native format). To calculate dn/dt, use the following:

Rich (BB code):
% Differentiates the vector n, i.e. dn/dt
ndiff = diff(n);
The same situation will apply for the variable p and dp/dt.

Can you draw up the code to initialise the simulation as above? We can look at the other stuff later.


Thread Starter

DJ Veni

Joined Jan 9, 2007
with knowing the Na and Nd concentrations we are able to find the majority carrier concentrations by the np=ni^2 relationship.

"you to get the n and p concentrations both minority and majority carrier concentrations with out the "influence of time" means, that the first thing I wish to do on matlab is to plot the concentration, the electric field and the band diagram on the y axis with the x axis denoting the length of the pn junction.
So time is not required in this first step.

Thread Starter

DJ Veni

Joined Jan 9, 2007
Attachment 1: N type semiconductor depletion region
p type semiconductor depletion region

Rich (BB code):
% Project Description
%Equilibrium Energy Band Diagram
%(Si,300K, nondegenerately doped step Junction)

NA=input('Enter p-side doping(cm^-3),NA=');
ND=input('Enter n-side doping(cm^-3),ND=');


Vbi=k*T*log((NA*ND)/ni^2);                                       % Built In Voltage 
xN=sqrt(2*KS*eO/q*NA*Vbi/(ND*(NA+ND)))                          % N-region depletion width   
xP=sqrt(2*KS*eO/q*ND*Vbi/(NA*(NA+ND)))                          % P-region depletion width
Vx1=(Vbi-q*ND.*(xN-x).^2/(2*KS*eO).*(x<=xN)).*(x>=0);            % Voltage at any point in n-region depletion region
Vx2=(.5*q*NA.*(xP+x).^2/(KS*eO).*(x>=-xP)).*(x<0);               % Voltage at any point in p-region depletion region 
VMAX=3;                                                          % Maximum Plot Voltage
EF=Vx(1)+VMAX/2-.026*reallog(NA/ni);                                  % Fermi Level

%Plot Diagram

axis([xleft xright 0 VMAX]);
hold on
plot([xleft xright],[EF EF],'w');
plot([0 0],[0.15 VMAX-0.5],'w--');
text(xleft*1.08,(-Vx(1)+EG/2+VMAX/2-.05),'Ec');                  %Plot label for Conduction Band as a Function of Vx and Doping at x<-xp
text(xright*1.02,(-Vx(200)+EG/2+VMAX/2-.05),'Ec');               %Plot label for Conduction Band as a Function of Vx and Doping at x>xn
text(xleft*1.08,(-Vx(1)-EG/2+VMAX/2-.05),'Ev');                  %Plot label for Valence Band as a Function of Vx and Doping at x<-xp
text(xright*1.02,(-Vx(200)-EG/2+VMAX/2-.05),'Ev');               %Plot label for Valence Band as a Function of Vx and Doping at x<-xp
text(xleft*1.08,(-Vx(1)+VMAX/2-.05),'Ei');                       %Plot label for Intrinsic Energy level.
hold off
put in Na as 4e15 and Nd as 3e14.

I need to apply depletion region approximation and then use the "!" the other stuff using differential equations.


Joined Nov 17, 2003

Further apologies for the lateness of my reply. The integrated differential function in Matlab will work by calculating the difference between adjacent elements in an array, which is derived with respect to a set of datum values, most commonly but not always time. So the key is to represent your differential function in terms of vectors - then it doesn't matter if you want to find a first, second, or fourteenth order derivative, the diff function will work. If you have multiple differentials then you can simply use the diff function as a stand alone on colaborate the solution in a new array:

e.g: d^2y/dt^2 + dx/dt = z

Can be deduced by:

Rich (BB code):
z = diff(y,2) + diff(x);
% Where x, y and z are all vectors derived with respect to t
Given I am not familiar with the "depletion region aproximation", I cannot comment specifically because I don't know how you intend to deduce the arrays upon which you are working. Sorry at this stage I cannot be more exact in my explanation (this is probably more to do with how long it is since I did semiconductor physics!)


Thread Starter

DJ Veni

Joined Jan 9, 2007

Well you just opened a door so that I can move foreward. Ill try it out and get back in a few days and let you know the results.

Dj Veni


Joined Nov 17, 2003
You could also have a look at the gradient function, however this too requires the input arguement to be a vector format. I have personally never used the gradient function, but from looking at its syntax in the Matlab help files, you may find it flexible enough to write you own diff function - certainly an option if you cannot get your data in a suitable format. I have also looked over at the Matlab file exchange to see if there are any other implementations of the diff function you can use, with little success.