This is version 0.5 There is a bug in this implementation, since hole concentration is constant in both side of the junction. Thanks to Urmimala for pointing this out.. In ver 1.0, we’ll provide complete documentation and will remove this bug.
The pn-junction is quite an interesting junction. If most the semiconductor theory can be tested around it, it also poses significant challenges in numerical method also. In this implementation in Scilab, we are solving Poisson equation to get the results when junction is in equilibrium. No part of this work used propriety softwares. Please wait for complete documentation of code.
Simulation Results
Click on image to enlarge it.
Scilab Functions
This is my top most function. Write a script which call this function and every other functions will be called in order.
[sourcecode language="matlab"]
// This function will discretized the domain and will create a grid points.
// Then we’ll discretise the equation there. Matrices
// will be created and saved in an binary file for further processing.
function [] = solveEquilibriumCondition()
// Open file to read the data.
// Check if file exists. If not run the script to create it.
ifExists = exists("./profilePN.dat");
if ifExists == 0 then
disp("File does not exist. Executing script to generate this file.");
exec("./profilePN.sci");
profilePN();
load("./profilePN.dat");
load("./myParameters.dat");
else
disp(‘File already exists. Loading …’);
load("./profilePN.dat");
load("./myParameters.dat");
disp("if there is any problem with this file. Delete it and re-execute this function/script.");
end;
delta_acc = 1E-5; // Preset the Tolerance
////////////////////////////////////////////////////////////////////////
/////////// EQUILIBRIUM SOLUTION PART BEGINS /////
////////////////////////////////////////////////////////////////////////
// Define the elements of the coefficient matrix for the internal nodes and
dx2 = delXdelX;
for i = 2: n_max-1
a(i) = 1/dx2;
c(i) = 1/dx2;
b(i) = -(2/dx2+exp(fi(i))+exp(-fi(i)));
f(i) = exp(fi(i)) – exp(-fi(i)) – dop(i) – fi(i)(exp(fi(i))+exp(-fi(i)));
end
// Define the elements of the coefficient matrix and initialize the forcing
// function at the ohmic contacts
a(1) = 0;
c(1) = 0;
b(1) = 1;
f(1) = fi(1);
a(n_max) = 0;
c(n_max) = 0;
b(n_max) = 1;
f(n_max) = fi(n_max);
// Start the iterative procedure for the solution of the linearized Poisson
// equation using LU decomposition method:
flag_conv = 0; // convergence of the Poisson loop
k_iter= 0;
while flag_conv == 0 then
k_iter = k_iter + 1;
alpha1(1) = b(1);
for i= 2 : n_max
beta1(i)=a(i)/alpha1(i-1);
alpha1(i)=b(i)-beta1(i)c(i-1);
end;
// Solution of Lv = f
v(1) = f(1);
for i = 2:n_max
v(i) = f(i) – beta1(i)v(i-1);
end;
// Solution of Ufi = v
temp = v(n_max)/alpha1(n_max);
delta(n_max) = temp – fi(n_max);
fi(n_max)=temp;
for i = (n_max-1):-1:1 //delta
temp = (v(i)-c(i)fi(i+1))/alpha1(i);
delta(i) = temp – fi(i);
fi(i) = temp;
end;
delta_max = 0;
for i = 1: n_max
xx = abs(delta(i));
if xx > delta_max then
delta_max=xx;
end;
//sprintf(‘delta_max = %d’,delta_max) //’k_iter = %d’,k_iter,’
end;
//delta_max=max(abs(delta));
// Test convergence and recalculate forcing function and
// central coefficient b if necessary
if delta_max < delta_acc then
flag_conv = 1;
else
for i = 2: n_max-1
b(i) = -(2/dx2 + exp(fi(i)) + exp(-fi(i)));
f(i) = exp(fi(i)) – exp(-fi(i)) – dop(i) – fi(i)(exp(fi(i)) + exp(-fi(i)));
end;
end;
end;
xx1(1) = delX1e4;
for i = 2:n_max-1
Ec(i) = dEc – Vtfi(i); //Values from the second Node%
ro(i) = -ni(exp(fi(i)) – exp(-fi(i)) – dop(i));
el_field1(i) = -(fi(i+1) – fi(i))Vt/(delXLdi);
el_field2(i) = -(fi(i+1) – fi(i-1))Vt/(2delXLdi);
n(i) = exp(fi(i));
p(i) = exp(-fi(i));
xx1(i) = xx1(i-1) + delXLdi1e4;
end;
Ec(1) = Ec(2);
Ec(n_max) = Ec(n_max-1);
xx1(n_max) = xx1(n_max-1) + delXLdi1e4;
el_field1(1) = el_field1(2);
el_field2(1) = el_field2(2);
el_field1(n_max) = el_field1(n_max-1);
el_field2(n_max) = el_field2(n_max-1);
nf = nni;
pf = pni;
ro(1) = ro(2);
ro(n_max) = ro(n_max-1);
// save all these values in a file.
save(‘equilibriumSol.dat’, Ec, Vt, fi, xx1, nf, pf, ro, n, p, el_field1, el_field2);
h = figure(1)
plot(xx1, Vtfi,’r’,’LineWidth’,2)
xlabel(‘x [um]’);
ylabel(‘Potential [eV]’);
title(‘Potential vs Position – at Equilibrium’);
xs2gif(1, ‘./figs/PotentialVsPos.gif’);
close(h);
h = figure(2)
plot(xx1, [el_field1, el_field2],’r’,’LineWidth’,2)
//plot(xx1, el_field2,’r’,’LineWidth’,2)
xlabel(‘x [um]’);
ylabel(‘Electric Field [V/cm]’);
title(‘Field Profile vs Position – at Equilibrium’);
xs2gif(2, ‘./figs/ElectricalFieldVsPos.gif’);
close(h);
h = figure(3)
plot2d1("onn", xx1, [nf, pf]);
xlabel(‘x [um]’);
ylabel(‘Electron & Hole Densities [1/cm^3]’);
title(‘Electron & Hole Densities vs Position – at Equilibrium’);
legend(‘n’,’p’);
xs2gif(3, ‘./figs/ElectronAndHolesDensity.gif’);
close(h);
h = figure(4)
plot(xx1, q_e*ro,’r’,’LineWidth’,2)
xlabel(‘x [um]’);
ylabel(‘Total Charge Density [C/cm^3]’);
title(‘Total Charge Density vs Position – at Equilibrium’);
xs2gif(4, ‘./figs/TotalChargeDensity.gif’);
close(h);
h = figure(5)
plot(xx1, Ec,’r’,’LineWidth’,2)
xlabel(‘x [um]’);
ylabel(‘Conduction Band Energy (eV)’);
title(‘Conduction Band vs Position – at Equilibrium’);
xs2gif(5, ‘./figs/ConductionBandDiagram.gif’);
close(h);
endfunction
[/sourcecode]
This function will profile the PN junction and return the doping and p and n distributions across the junction
[sourcecode language="matlab"]
// This function will profile the PN junction and return the doping and p and n
// distributions across the junction.
// Open file to read the data.
// Check if file exists. If not run the script to create it.
function [] = profilePN()
ifExists = exists("./myParameters.dat");
if ifExists == 0 then
disp("File does not exist. Executing script to generate this file.");
exec("./calcParameters.sci");
calcParameters();
load("./myParameters.dat");
else
disp('File already exists. Loading ...');
load("./myParameters.dat");
disp("if there is any problem with this file. Delete it and re-execute this function/script.");
end;
// Now we need to discretise the boundries.
// Set grid size based on the extrinsic Debye Lengths.
if Ldn < Ldp then
delX = Ldn/10;
else
delX = Ldp/10;
end;
// Length of the pn junction.
if Wp > Wn then
length_pn = 30*Wn;
else
length_pn = 30*Wp;
end;
// We have grid size, We have length. Calculate the total points.
n_max = floor(length_pn/delX);
delX = delX/Ldi; // renormalize lengths with Ldi.
// Set up doping C(x) = Nd(x) - Na(x) that is normalised with ni
// Half of the point are p type, rests are n-types.
for i = 1:n_max
if i <= n_max/2 then
dop(i) = -Na/ni; // p type
elseif i > n_max/2 then
dop(i) = Nd/ni; // n type
end
end
// Initialize the potential based on the requirement of charge neutrality
// throughout the whole structure.
for i = 1: n_max
tempZ = 0.5*dop(i);
if tempZ > 0 then
tempX = tempZ * (1 + sqrt(1 + 1/(tempZ*tempZ)));
elseif tempZ < 0 then
tempX = tempZ * (1 - sqrt( 1 + 1/(tempZ*tempZ)));
end;
fi(i) = log(tempX);
n(i) = tempX;
p(i) = 1/tempX;
end;
save('profilePN.dat', fi, p, n, dop, length_pn, n_max, delX);
endfunction
[/sourcecode]
Here is the part which calculates the parameters to be used by above posted functions. [sourcecode language="matlab"]
// This function will calculate the parameters. We'll store these parameters
// in a file. We'll use this file to reuse the parametes for further processing.
// (c) Dilawar, 2010
// dilawar@ee.iitb.ac.in
// www.ee.iitb.ac.in/student/~dilawar
// GNU - GPL
// This function should return a ascii file containing parameters which are to
// be used later.
function [] = calcParameters()
T = 300; // Substrate is Silicon. Temp 300 K
u_e = 1350; // cm*cm/V/s Mobility of electron in Silicon
u_p = 480; // cm*cm/V/s Mobility of holes in Silicon
D_n = 35; // cm*cm/s Diffusion coeff.
D_p = 12.4; // cm*cm/s Diffusion coeff.
q_e = 1.6e-19; // Charge on an electron.
kt_q = 0.0259;
kb = 1.38e-23; // Boltzmann Constant JK^-1
eps = 1.05e-12; // Permittivity of undopes silicon.
ni = 1.5e10; // Intrinsic carrier concentration
Vt = kb*T/q_e; // Threshold voltages.
RNc = 2.8e20;
dEc = Vt*log(RNc/ni);
tauN0 = 0.1e-6;
tauP0 = 0.1e-6;
mu_n0 = 1500;
mu_p0 = 1000;
// Define doping values.
Na = 1e16;
Nd = 1e18;
// Calculate parameters.
Vbi = Vt*log(Na*Nd/(ni*ni)); // Built-in voltage.
W = sqrt(2*eps*(Na + Nd)*Vbi/(q_e*Na*Nd) ); // Depletion width
Wn = W*sqrt(Na/(Na + Nd)); // Depletion width at n side.
Wp = W*sqrt(Nd/(Na + Nd)); // Depletion width at p side
Wone = sqrt(2*eps*Vbi/(q_e*Na)); //
E_p = q_e*Nd*Wn/eps;
Ldn = sqrt(eps*Vt/(q_e*Nd));
Ldp = sqrt(eps*Vt/(q_e*Na));
Ldi = sqrt(eps*Vt/(q_e*ni));
// We'd like to save our paramters in some file.
save('myParameters.dat', q_e, Na, Nd, Vbi, W, Wn, Wp, E_p, Ldn, Ldp, Ldi, ni, dEc, Vt, kt_q, mu_n0, mu_p0);
endfunction;
[/sourcecode]
// This function will discretize the domain and will create a grid points. // Then we'll discretise the equation there. Matrices // will be created and saved in an binary file for further processing. function [] = solveEquilibriumCondition() // Open file to read the data. // Check if file exists. If not run the script to create it. ifExists = exists("./profilePN.dat"); if ifExists == 0 then disp("File does not exist. Executing script to generate this file."); exec("./profilePN.sci"); profilePN(); load("./profilePN.dat"); load("./myParameters.dat"); else disp('File already exists. Loading ...'); load("./profilePN.dat"); load("./myParameters.dat"); disp("if there is any problem with this file. Delete it and re-execute this function/script."); end; delta_acc = 1E-5; // Preset the Tolerance //////////////////////////////////////////////////////////////////////// /////////// EQUILIBRIUM SOLUTION PART BEGINS ///// //////////////////////////////////////////////////////////////////////// // Define the elements of the coefficient matrix for the internal nodes and dx2 = delX*delX; for i = 2: n_max-1 a(i) = 1/dx2; c(i) = 1/dx2; b(i) = -(2/dx2+exp(fi(i))+exp(-fi(i))); f(i) = exp(fi(i)) - exp(-fi(i)) - dop(i) - fi(i)*(exp(fi(i))+exp(-fi(i))); end // Define the elements of the coefficient matrix and initialize the forcing // function at the ohmic contacts a(1) = 0; c(1) = 0; b(1) = 1; f(1) = fi(1); a(n_max) = 0; c(n_max) = 0; b(n_max) = 1; f(n_max) = fi(n_max); // Start the iterative procedure for the solution of the linearized Poisson // equation using LU decomposition method: flag_conv = 0; // convergence of the Poisson loop k_iter= 0; while flag_conv == 0 then k_iter = k_iter + 1; alpha1(1) = b(1); for i= 2 : n_max beta1(i)=a(i)/alpha1(i-1); alpha1(i)=b(i)-beta1(i)*c(i-1); end; // Solution of Lv = f v(1) = f(1); for i = 2:n_max v(i) = f(i) - beta1(i)*v(i-1); end; // Solution of U*fi = v temp = v(n_max)/alpha1(n_max); delta(n_max) = temp - fi(n_max); fi(n_max)=temp; for i = (n_max-1):-1:1 //delta temp = (v(i)-c(i)*fi(i+1))/alpha1(i); delta(i) = temp - fi(i); fi(i) = temp; end; delta_max = 0; for i = 1: n_max xx = abs(delta(i)); if xx > delta_max then delta_max=xx; end; //sprintf('delta_max = %d',delta_max) //'k_iter = %d',k_iter,' end; //delta_max=max(abs(delta)); // Test convergence and recalculate forcing function and // central coefficient b if necessary if delta_max < delta_acc then flag_conv = 1; else for i = 2: n_max-1 b(i) = -(2/dx2 + exp(fi(i)) + exp(-fi(i))); f(i) = exp(fi(i)) - exp(-fi(i)) - dop(i) - fi(i)*(exp(fi(i)) + exp(-fi(i))); end; end; end; xx1(1) = delX*1e4; for i = 2:n_max-1 Ec(i) = dEc - Vt*fi(i); //Values from the second Node% ro(i) = -ni*(exp(fi(i)) - exp(-fi(i)) - dop(i)); el_field1(i) = -(fi(i+1) - fi(i))*Vt/(delX*Ldi); el_field2(i) = -(fi(i+1) - fi(i-1))*Vt/(2*delX*Ldi); n(i) = exp(fi(i)); p(i) = exp(-fi(i)); xx1(i) = xx1(i-1) + delX*Ldi*1e4; end; Ec(1) = Ec(2); Ec(n_max) = Ec(n_max-1); xx1(n_max) = xx1(n_max-1) + delX*Ldi*1e4; el_field1(1) = el_field1(2); el_field2(1) = el_field2(2); el_field1(n_max) = el_field1(n_max-1); el_field2(n_max) = el_field2(n_max-1); nf = n*ni; pf = p*ni; ro(1) = ro(2); ro(n_max) = ro(n_max-1); // save all these values in a file. save('equilibriumSol.dat', Ec, Vt, fi, xx1, nf, pf, ro, n, p, el_field1, el_field2); h = figure(1) plot(xx1, Vt*fi,'r','LineWidth',2) xlabel('x [um]'); ylabel('Potential [eV]'); title('Potential vs Position - at Equilibrium'); xs2gif(1, './figs/PotentialVsPos.gif'); close(h); h = figure(2) plot(xx1, [el_field1, el_field2],'r','LineWidth',2) //plot(xx1, el_field2,'r','LineWidth',2) xlabel('x [um]'); ylabel('Electric Field [V/cm]'); title('Field Profile vs Position - at Equilibrium'); xs2gif(2, './figs/ElectricalFieldVsPos.gif'); close(h); h = figure(3) plot2d1("onn", xx1, [nf, pf]); xlabel('x [um]'); ylabel('Electron & Hole Densities [1/cm^3]'); title('Electron & Hole Densities vs Position - at Equilibrium'); legend('n','p'); xs2gif(3, './figs/ElectronAndHolesDensity.gif'); close(h); h = figure(4) plot(xx1, q_e*ro,'r','LineWidth',2) xlabel('x [um]'); ylabel('Total Charge Density [C/cm^3]'); title('Total Charge Density vs Position - at Equilibrium'); xs2gif(4, './figs/TotalChargeDensity.gif'); close(h); h = figure(5) plot(xx1, Ec,'r','LineWidth',2) xlabel('x [um]'); ylabel('Conduction Band Energy (eV)'); title('Conduction Band vs Position - at Equilibrium'); xs2gif(5, './figs/ConductionBandDiagram.gif'); close(h); endfunction