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
% Version : 1.1
% 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
% Category : live script
% To Do : -
% -------------------------------------------------------------------------
% 1:01: initial version [PG]
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
Table of Contents
Poisson Equation for Electrostatics Newton-Raphson Method Finite Difference Formulation Discretization of the 1D Poisson Equation Declaration of Constants Grid Specifications Doping Profile Newton-Raphson Parameters Implementation of a Newton-Raphson Assisted Finite Difference 1D Poisson Solver Plotting the Results Doping Profile Charge Carrier Concentrations Charge Density Internal Electric Field Build-in Potential Band Diagram References

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)
with , , , and
. (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
% (temperature-dependent band-gap and intrinsic carrier concentration is not implemented)
T = 300; % in degree K
% permittivity of vacuum
eps_0 = 8.854187813e-12; % in As/Vm
% relative permittivity
eps_r = 12;
% absolute permittivity
epsilon = eps_0 * eps_r; % in As/Vm
% Boltzmann constant
kB = 1.380658e-23; % in J/K
% elementary charge
q = 1.60217733e-19; % in As
% intrinsic carrier concentration of SI
ni = 1.45e16; % in 1/m^3
% band gap energy of SI
Eg = 1.12; % in eV

Grid Specifications

% lower grid bound
xmin = str2double("-4e-6"); % in m
% upper grid bound
xmax = str2double("4e-6"); % in m
% number of grid points
N = round(str2double("1000"));

Doping Profile

doping_profile = "abrupt_pnp";

Newton-Raphson Parameters

maxIter = round(str2double("1000"));
alpha = 1e-6;

Implementation of a Newton-Raphson Assisted Finite Difference 1D Poisson Solver

% grid points (sampling)
x = linspace(xmin,xmax,N).';
% step-size
dx = x(2)-x(1);
% grid extension (needed for central differences)
xg = [ min(x)-dx;x;max(x)+dx ];
% load doping profile
[NA,ND] = doping_profiles(doping_profile,x);
% initializing electron concentration
n = ND;
% initializing hole concentration
p = NA;
% initializing electric potential
V = zeros(size(xg));
Vinit = band_initialization(NA,ND,ni,q,kB,T);
V(2:N+1) = Vinit;
V(1) = V(2); % boundary
V(end) = V(end-1); % boundary
V_old = V;
% figure;
% plot(x,V(2:N+1))
% loop for Newton-Raphson Iterations
for k=1:maxIter
% updating charge carrier concentrations
n = ni*exp(q*V(2:N+1)/(kB*T));
p = ni*exp(-q*V(2:N+1)/(kB*T));
% charge density
% (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
delta_V = M_k\R_k;
% 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;
% stopping criteria
if( norm( (V - V_old), "inf" ) < alpha )
break;
end
V_old = V;
end
% calculating the electric field
E = -diff(V(1:N+1))/dx;
% build in voltage
Vbi = abs(V(end)-V(1));
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

set(groot, ...
'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')
xlabel('$x$ in $m$')
ylabel('$N(x)$ in $1$/$m^3$')
legend_obj = legend();
legend_obj.Location = 'northeastoutside';
hold off;

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')
xlabel('$x$ in $m$')
ylabel('$n(x),p(x)$ in $1$/$m^3$')
legend_obj = legend();
legend_obj.Location = 'northeastoutside';
hold off;

Charge Density

figure, hold on, grid minor;
plot(x,rho,'b','DisplayName','$\rho$','LineWidth',1.5)
title('charge density')
xlabel('$x$ in $m$')
ylabel('$\rho(x)$ in $As$/$m^3$')
legend_obj = legend();
legend_obj.Location = 'northeastoutside';
hold off;

Internal Electric Field

figure, hold on, grid minor;
plot(x,E,'b','DisplayName','$E$','LineWidth',1.5)
title('electric field')
xlabel('$x$ in $m$')
ylabel('$E(x)$ in $V$/$m$')
legend_obj = legend();
legend_obj.Location = 'northeastoutside';
hold off;

Build-in Potential

figure, hold on, grid minor;
plot(x,V(2:N+1),'b','DisplayName','$V$','LineWidth',1.5)
title('potential')
xlabel('$x$ in $m$')
ylabel('$V(x)$ in $V$')
legend_obj = legend();
legend_obj.Location = 'northeastoutside';
hold off;

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)
title('band diagram')
xlabel('$x$ in $m$')
ylabel('energy in $eV$')
legend_obj = legend();
legend_obj.Location = 'northeastoutside';
hold off;
function V = band_initialization(NA,ND,ni,q,kB,T)
% if needed, add bias voltage to this function
dN = ND - NA;
V = zeros(size(dN));
for i=1:length(V)
if( abs(dN) >= ni )
V(i) = sign(dN(i))*kB*T/q*log(abs(dN(i))/ni);
end
end
end
function [NA,ND] = doping_profiles(str,x)
switch str
case 'abrupt'
NA = 1E14*ones(size(x)); % 1/cm^3
NA = NA*1e6; % 1/m^3
ND = [ zeros(size(x(x<0))); ones(size(x(x>=0))) ] * 2E14; % 1/cm^3
ND = ND*1e6; % 1/m^3
case 'linear'
alpha = 1e19*100e6; % 1/m^4
NA = zeros(size(x));
ND = zeros(size(x));
for i=1:length(x)
if( x(i) <= 0 )
NA(i) = -alpha*x(i);
elseif( x(i) > 0 )
ND(i) = alpha*x(i);
end
end
case 'hyper_abrupt'
NA=1E15*exp(0.5*(x/1e-6)).*heaviside(-x); % 1/cm^3
NA = NA*1e6; % 1/m^3
ND=1E15*exp(0.5*(-x/1e-6)).*heaviside(x); % 1/cm^3
ND = ND*1e6; % 1/m^3
case 'abrupt_pnp'
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
case 'abrupt_npn'
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
case 'abrupt_pnpn'
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
case 'abrupt_npnp'
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
case 'arbitrary'
[NA,ND] = myProfile(x);
otherwise
error(['[ERROR]: ',str,' is not a valid doping profile!'])
end
end
function [NA,ND] = myProfile(x)
NA = 1E14*ones(size(x)); % 1/cm^3
NA = NA*1e6; % 1/m^3
ND = exp( -(x/1e-6).^2 )*1e15; % 1/cm^3
ND = ND*1e6; % 1/m^3
end

References