Semiconductor Physics P-N junction

Discussion in 'Programmer's Corner' started by DJ Veni, Jan 11, 2007.

  1. DJ Veni

    Thread Starter New Member

    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.

  2. Dave

    Retired Moderator

    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.

  3. DJ Veni

    Thread Starter New Member

    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.
    • Fig.jpg
      File size:
      74.6 KB
  4. Dave

    Retired Moderator

    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:

    Code ( (Unknown Language)):
    1. 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:

    Code ( (Unknown Language)):
    1. 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:

    Code ( (Unknown Language)):
    1. plot(x,e)
    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.

  5. DJ Veni

    Thread Starter New Member

    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 : [​IMG]

    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.
  6. Dave

    Retired Moderator

    Nov 17, 2003
    Please check your Private Messages.

  7. Dave

    Retired Moderator

    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:

    Code ( (Unknown Language)):
    1. % Lines preceded with one of these => is a comment
    2. % Set ni equal to 1x10^10
    3. 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:

    Code ( (Unknown Language)):
    1. 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:

    Code ( (Unknown Language)):
    1. % Define t from 0 to 60 second at intervals of 0.5 seconds
    2. 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:

    Code ( (Unknown Language)):
    1. % Differentiates the vector n, i.e. dn/dt
    2. 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.

  8. DJ Veni

    Thread Starter New Member

    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.
  9. DJ Veni

    Thread Starter New Member

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

    Code ( (Unknown Language)):
    1. % Project Description
    2. %Equilibrium Energy Band Diagram
    3. %(Si,300K, nondegenerately doped step Junction)
    5.                                                      %Constants
    6. q=1.6e-19;
    7. KS=11.8;
    8. ni=1.0e10;
    9. EG=1.12;
    10. T=300;
    11. k=8.617e-5;
    12. eO=8.854e-14;
    13. xleft=-5e-4;
    14. xright=-xleft;
    15. NA=input('Enter p-side doping(cm^-3),NA=');
    16. ND=input('Enter n-side doping(cm^-3),ND=');
    18.                                                      %Computations
    20. Vbi=k*T*log((NA*ND)/ni^2);                                       % Built In Voltage
    21. xN=sqrt(2*KS*eO/q*NA*Vbi/(ND*(NA+ND)))                          % N-region depletion width  
    22. xP=sqrt(2*KS*eO/q*ND*Vbi/(NA*(NA+ND)))                          % P-region depletion width
    23. x=linspace(xleft,xright,200);                                      
    24. Vx1=(Vbi-q*ND.*(xN-x).^2/(2*KS*eO).*(x<=xN)).*(x>=0);            % Voltage at any point in n-region depletion region
    25. Vx2=(.5*q*NA.*(xP+x).^2/(KS*eO).*(x>=-xP)).*(x<0);               % Voltage at any point in p-region depletion region
    26. Vx=Vx1+Vx2;                                                      
    27. VMAX=3;                                                          % Maximum Plot Voltage
    28. EF=Vx(1)+VMAX/2-.026*reallog(NA/ni);                                  % Fermi Level
    30. %Plot Diagram
    32. close
    33. plot(x,-Vx+EG/2+VMAX/2);
    34. axis([xleft xright 0 VMAX]);
    35. axis('off');
    36. hold on
    37. plot(x,-Vx-EG/2+VMAX/2);                                          
    38. plot(x,-Vx+VMAX/2,'w');
    39. plot([xleft xright],[EF EF],'w');
    40. plot([0 0],[0.15 VMAX-0.5],'w--');
    41. %plot(xP,EF,'XP');
    42. %plot(-xN,EF,'XN');
    43. 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
    44. 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
    45. 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
    46. 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
    47. text(xleft*1.08,(-Vx(1)+VMAX/2-.05),'Ei');                       %Plot label for Intrinsic Energy level.
    48. text(xright*1.02,EF-.05,'EF');
    49. set(gca,'DefaultTextUnits','Normalized')
    50. text(.18,0,'p-side');
    51. text(.47,0,'x=0');
    52. text(.75,0,'n-side');
    53. set(gca,'defaultTextUnits','Data')
    54. 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.
  10. Dave

    Retired Moderator

    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:

    Code ( (Unknown Language)):
    1. z = diff(y,2) + diff(x);
    2. % 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!)

  11. DJ Veni

    Thread Starter New Member

    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
  12. Dave

    Retired Moderator

    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.