Finite Difference based 1D Nonlinerar Possion Equation Solver for pn-Junctions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Name : pn_poisson_calc.mlx
% Author : Patrick Gsoels (11771814)
% Contact : patrick.gsoels@student.tugraz.at
% Creation Date : 16.07.2022
% Description : This live script computes the charge carrier
% concentrations, charge-density, electric field,
% elctric potential and band-structure of arbitrary
% pn-junction doping profiles.
% Language : MATLAB R2020b
% -------------------------------------------------------------------------
% 1:01: initial version [PG]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Poisson Equation for Electrostatics
In electrodynamics Gauss's law states the relation between the electric displacement field and the free charge density ρ within a volume. The differential form of Gauss's law can be written as . (1) For a homogenious medium the following relation between electric displacement field and electric field holds: (2) In electrostatics we do not have any time-dependent fields. Hence, Faraday's law of induction will simplify to the following expression:
(3) The simplified version of Faraday's law of induction states that the electric field is free of curl. Thus, the electric field is a conservative field for which one can define a potential. (4) Inserting Equation (2) and (4) into Equation (1) will give us the Poisson equation for electrostatics.
(5) For a constant dielectric media Equation (5) simplifies to
. (6) pn-junctions do typically have doping profiles which intentionally vary along the latteral axis which is perpendicular to the junction interface. Thus, the Poisson equation from (6) reduces to one spacial coordinate.
. (7) Due to the interconnection of differently doped materials, charge carriers diffuse accros the pn-junction. At the same time an internal electric field establishes due to the movement of charge carriers. These facts lead to the necessity of introducing a potential-dependent charge density ρ, which reformulates the Poisson equation to
. (8) According to [1], assuming the dopants to be totally ionised, the charge density inside a semiconductor is given by
. (9) According to [2], the electron and hole density can be formulated as (10) and
. (11) Inserting Equation (10) and (11) into (9) will provide us with the full formulation of the potential-dependent charge density ρ.
(12) with being the intrinsic carrier concentration of the semiconductor, which is defined as . (13) Bear in mind that the formulars for and are derived by using the Boltzmann approximation instead of the real fermi-dirac distribution. Hence, the condition has to hold. With the help of some trigonometric identity (), Equation (12) can be reformulated as . (14) Since our one-dimensional Poisson equation is non-linear, one will need at some point the partial derivative of the charge density ρ with respect to the electric potential V.
(15) Newton-Raphson Method
The Newton-Raphson method is mathematically the most preferred of the several known methods for the solution of systems of nonlinear equations. It is a successive approximation technique based on an initial estimate of the solution and the Taylor series expansion [1]. Consider the set of equations
(16) A Taylor approximation of (16), which is terminated after the linear term, leads to
. (17) . (18) with
(19) being the Jaccobian.
After rearranging (17), one gets the expression for the Newton-Raphson step
(20) and the ()-th iterate . (21) Finite Difference Formulation
The finite difference formulation can be again derived from the Taylor approximation of as depicted below. (22) [forward difference] (23) (24) [backward difference] (25) One gets the central difference formulation for the first order derivative if one subtracts Equation (24) from (22) like shown below.
(26) [central difference] (27) The second order derivative can be approximated by inserting Equation (23) and (25) into the expression depicted below.
(28) Discretization of the 1D Poisson Equation
With the aid of the finite difference formulation one can approximate Equation (8), which results in
, . (29) After rearranging Equation (29) one gets the non-linear algebraic difference formulation of the Poisson Equation.
, . (30) By using the Newton-Raphson technique for solving non-linear systems of equations, the non-linear algebraic difference equation (30) formulates to
. (31) The known parts of Equation (31) will define our right hand-side and the unknown parts will define the left-hand side of the system of equations which needs to be solved for . (32) By defining the right-hand side of Equation (32) as
(33) and the left-hand side as
. (34) With Equation (33) and (34) one can easily write Equation (32) in matrix form
(35) . (36) Finally, one has to solve Equation (35) for the Newton-Raphson step in order to generate the next iterate of the potential as defined by Equation (21). (37) Declaration of Constants
% (temperature-dependent band-gap and intrinsic carrier concentration is not implemented)
eps_0 = 8.854187813e-12; % in As/Vm
epsilon = eps_0 * eps_r; % in As/Vm
kB = 1.380658e-23; % in J/K
q = 1.60217733e-19; % in As
% intrinsic carrier concentration of SI
Grid Specifications
xmin = str2double("-4e-6"); % in m
xmax = str2double("4e-6"); % in m
N = round(str2double("1000"));
Doping Profile
doping_profile = "abrupt_pnp";
Newton-Raphson Parameters
maxIter = round(str2double("1000"));
Implementation of a Newton-Raphson Assisted Finite Difference 1D Poisson Solver
x = linspace(xmin,xmax,N).';
% grid extension (needed for central differences)
xg = [ min(x)-dx;x;max(x)+dx ];
[NA,ND] = doping_profiles(doping_profile,x);
% initializing electron concentration
% initializing hole concentration
% initializing electric potential
Vinit = band_initialization(NA,ND,ni,q,kB,T);
V(end) = V(end-1); % boundary
% loop for Newton-Raphson Iterations
% updating charge carrier concentrations
n = ni*exp(q*V(2:N+1)/(kB*T));
p = ni*exp(-q*V(2:N+1)/(kB*T));
% (assuming that the dopants are totally ionized)
rho = q*(p - n + ND - NA);
% charge density change due to the internal potential
drho_dV = -2*q^2*ni/(kB*T) * cosh(q*V(2:N+1)/(kB*T));
% finite difference formulation for the potential
d2V_dx2 = ( V(1:N,1) - 2*V(2:N+1,1) + V(3:N+2,1) )/(dx^2);
% non-linear equation system
Mjj_k = 2/(dx^2) - 1/epsilon*drho_dV;
M_off = 1/(dx^2)*ones(N,1);
M_k = spdiags([-M_off Mjj_k -M_off],-1:1,N,N);
R_k = 1/epsilon*rho + d2V_dx2;
% solving for potential changes
% Newton-Raphson Step (updating the potential)
%V(2:N+1) = V(2:N+1)+(0.98^k)*delta_V;
V(2:N+1) = V(2:N+1)+delta_V;
if( norm( (V - V_old), "inf" ) < alpha )
% calculating the electric field
disp(['Number of Newton-Raphson Iterations: ',num2str(k)])
Number of Newton-Raphson Iterations: 6
disp(['build in voltage: Vbi = ',num2str(Vbi),'V'])
build in voltage: Vbi = 0.041607V
Plotting the Results
'DefaultTextInterpreter', 'LaTeX',...
'DefaultAxesTickLabelInterpreter', 'LaTeX',...
'DefaultLegendInterpreter', 'LaTeX',...
'DefaultAxesFontName', 'LaTeX',...
'DefaultFigureColor', 'w', ...
'DefaultAxesLineWidth', 0.5, ...
'DefaultAxesXColor', 'k', ...
'DefaultAxesYColor', 'k', ...
'DefaultAxesFontUnits', 'points', ...
'DefaultAxesFontSize', 15, ...
'DefaultAxesFontName', 'Helvetica', ...
'DefaultLineLineWidth', 1.5, ...
'DefaultTextFontUnits', 'Points', ...
'DefaultTextFontSize', 20, ...
'DefaultTextFontName', 'Helvetica', ...
'DefaultAxesBox', 'on', ...
'DefaultAxesTickLength', [0.01 0.01],...
'defaultAxesXGrid','on',...
'defaultAxesYGrid','on',...
'defaultfigureposition',[50 50 1080 400],...
'defaultAxesTitleFontSize',1.3);
Doping Profile
% plot chosen doping profile
figure, hold on, grid minor;
plot(x,NA,'b','DisplayName','$N_A$','LineWidth',1.5)
plot(x,ND,'r','DisplayName','$N_D$','LineWidth',1.5)
plot(x,(ND-NA),'g--','DisplayName','$N_D - N_A$','LineWidth',1.5)
title('doping concentrations')
ylabel('$N(x)$ in $1$/$m^3$')
legend_obj.Location = 'northeastoutside';
Charge Carrier Concentrations
figure, hold on, grid minor;
plot(x,p,'b','DisplayName','$p$','LineWidth',1.5)
plot(x,n,'r','DisplayName','$n$','LineWidth',1.5)
title('carrier concentrations')
ylabel('$n(x),p(x)$ in $1$/$m^3$')
legend_obj.Location = 'northeastoutside';
Charge Density
figure, hold on, grid minor;
plot(x,rho,'b','DisplayName','$\rho$','LineWidth',1.5)
ylabel('$\rho(x)$ in $As$/$m^3$')
legend_obj.Location = 'northeastoutside';
Internal Electric Field
figure, hold on, grid minor;
plot(x,E,'b','DisplayName','$E$','LineWidth',1.5)
ylabel('$E(x)$ in $V$/$m$')
legend_obj.Location = 'northeastoutside';
Build-in Potential
figure, hold on, grid minor;
plot(x,V(2:N+1),'b','DisplayName','$V$','LineWidth',1.5)
legend_obj.Location = 'northeastoutside';
Band Diagram
figure, hold on, grid minor;
plot(x,-V(2:N+1)+Eg/2,'r','DisplayName','conduction band','LineWidth',1.5)
plot(x,-V(2:N+1)-Eg/2,'b','DisplayName','valence band','LineWidth',1.5)
legend_obj.Location = 'northeastoutside';
function V = band_initialization(NA,ND,ni,q,kB,T)
% if needed, add bias voltage to this function
V(i) = sign(dN(i))*kB*T/q*log(abs(dN(i))/ni);
function [NA,ND] = doping_profiles(str,x)
NA = 1E14*ones(size(x)); % 1/cm^3
ND = [ zeros(size(x(x<0))); ones(size(x(x>=0))) ] * 2E14; % 1/cm^3
alpha = 1e19*100e6; % 1/m^4
NA=1E15*exp(0.5*(x/1e-6)).*heaviside(-x); % 1/cm^3
ND=1E15*exp(0.5*(-x/1e-6)).*heaviside(x); % 1/cm^3
NA = ( heaviside(-x-1e-6)*5E15 + heaviside(x+1e-6)*1E15 )*1e6; % 1/m^3
ND = ( heaviside(x+1e-6)*3E15 - heaviside(x-1e-6)*3E15 )*1e6; % 1/m^3
ND = ( heaviside(-x-1e-6)*5E15 + heaviside(x+1e-6)*1E15 )*1e6; % 1/m^3
NA = ( heaviside(x+1e-6)*3E15 - heaviside(x-1e-6)*3E15 )*1e6; % 1/m^3
NA = ( heaviside(-x-2e-6)*1E18 + ( heaviside(x-2e-6) - heaviside(x-3e-6) )*0.8E18 )*1e6; % 1/m^3
ND = ( ( heaviside(x+2e-6) - heaviside(x-2e-6) )*5E17 + heaviside(x-3e-6)*1.2E18 )*1e6; % 1/m^3
ND = ( heaviside(-x-2e-6)*1E18 + ( heaviside(x-2e-6) - heaviside(x-3e-6) )*0.8E18 )*1e6; % 1/m^3
NA = ( ( heaviside(x+2e-6) - heaviside(x-2e-6) )*5E17 + heaviside(x-3e-6)*1.2E18 )*1e6; % 1/m^3
error(['[ERROR]: ',str,' is not a valid doping profile!'])
function [NA,ND] = myProfile(x)
NA = 1E14*ones(size(x)); % 1/cm^3
ND = exp( -(x/1e-6).^2 )*1e15; % 1/cm^3
References
- [1] Jabr, R.A. & Hamad, Mohammed & Mohanna, Yasser. (2007). Newton-Raphson Solution of Poisson's Equation in a Pn Diode. International Journal of Electrical Engineering Education. 44. 23-33. 10.7227/IJEEE.44.1.3. https://www.researchgate.net/publication/233646108_Newton-Raphson_Solution_of_Poisson%27s_Equation_in_a_Pn_Diode
- [2] SZE, S. M. (1985). Semiconductor devices, physics and technology. New York, Wiley.
- [3] F. R. Shapiro, "The numerical solution of Poisson's equation in a pn diode using a spreadsheet," in IEEE Transactions on Education, vol. 38, no. 4, pp. 380-384, Nov. 1995, doi: 10.1109/13.473161.
- [4] Chapra, Steven C, and Raymond P. Canale. Numerical Methods for Engineers. Boston: McGraw-Hill Higher Education, 2006. Print.
- [5] Lecture 9 (Computational Electromagnetics) - Finite Difference Method, Dr. Raymond Rumpf, https://www.youtube.com/watch?v=v-exTNOSG3g