% quasi-2D Poisson-Continuity Equation Solver
% Author: Xiao Wang
% The University of Texas at Austin,
% J.J. Pickle Research Center,
% Bldg. #160, 10100 Burnet Rd.
% Austin, TX 78758
% Copyright (c) 2022, Xiao Wang (please also see the license)
% common constant in use
e = 1.6E-19;
h = 6.63E-34;
hbar = h/2/pi;
kb = 1.38E-23;
E0 = 8.85*10^(-12);
m0 = 9.1E-31;
% semiconductor material parameters
% (a-IGZO as an example)
Ec = 0;
me = 0.34*m0;
w0 = 20E-3*e/hbar;
a = 0.33E-9;
Beta = 2;
gs = 2;
gv = 1;
Es = 10;
Esinf = 4;
Eox = 3.9;
%trap profile
Nt = 5*10^16;
gamma = 1;
%device profile
%gate capacitance
tox = 100E-9;
Cox = E0*Eox/tox;
W = 10E-6;
%channel length
L = 1E-6;
nL = 50;
Lpo = linspace(0,L,nL)’;
dL = Lpo(2)-Lpo(1);
%operating conditions
Id = 1E-5;
Vd = zeros(length(Id),1);
Vov = 20*ones(length(Vd),1);
Vov0 = Vov(1);
T = 300;
%calculate channel thickness (quantum well width)
Ft0 = Cox.*Vov0/(E0*Es);
an1 = -2.338;
E1 = -(e^2*Ft0.^2*hbar^2/(2*me)).^(1/3)*an1;
zav = 2*E1./(3*e*Ft0);
d = 2*zav;
%% solve Vd for certain Id
for iVd = 1:length(Id)
%% mobility and electric field at source
%Pmtr at source end
syms En Efs
g_band = gs*gv*me/(2*pi*hbar^2);
DOS_band = g_band.*heaviside(En-Ec);
Tta = Nt/(g_band)/kb;
g_trap = Nt/(kb*Tta)*exp(En/(kb*Tta));
DOS_trap = g_trap.*heaviside(Ec-En);
DOS_tot = DOS_band + DOS_trap;
FE = 1./(1+exp((En-Efs)/(kb*T)));
ntotint = FE*DOS_tot;
ntot0 = Cox*(Vov0)/e;
eqn = ntot0 == vpaintegral(ntotint,En,-2*e,2*e);
Efsol0 = double(vpasolve(eqn,Efs));
Eint = e*linspace(-5,5,1000);
FEint = 1./(1+exp((Eint-Efsol0)/(kb*T)));
ntrapint = FEint.*Nt/(kb*Tta).*exp(Eint/(kb*Tta)).*heaviside(Ec-Eint);
nbandint = FEint.*gs*gv*me/(2*pi*hbar^2).*heaviside(Eint-Ec);
ntrap0 = trapz(Eint,ntrapint);
nband0 = ntot0-ntrap0;
Pmtr0 = 1-ntrap0./ntot0;
%band mobility at source
%trapped carrier scattering at source end
kav = sqrt(2*kb*T*me)/hbar; %non-degenerate
kf = sqrt(2*pi*nband0); %degenerate
theta = linspace(0,pi)’;
S1 = e^2*nband0/(2*kb*T*E0*Es); %non-degenerate
S2 = 2*e^2*me/(4*pi*E0*Es*hbar^2); %degenerate
S0 = S1.*heaviside(S2-S1)+S2.*heaviside(S1-S2);
keff = kav.*heaviside(S2-S1)+kf.*heaviside(S1-S2);
Ibbint = (sin(theta)).^2./(sin(theta)+S0/keff).^2;
Ibb = trapz(theta,Ibbint);
tao_tc0 = 1./(e^4*(ntrap0)/d*me/(8*pi*(E0*Es)^2*hbar^3*(keff).^3)*Ibb);
%PO_phonon scattering at source end
Ediff = 1/(1/Esinf-1/Es);
Np = 1/(exp(hbar*w0/(kb*T))-1);
tao_po0 = 4*pi*E0*Ediff*hbar^2/(e^2*w0*me*d)/Np;
%band mobility at source end
taoband0 = 1/(1/tao_po0+1/tao_tc0);
miuband0 = taoband0*e/me;
%TRF factor at source end
vth0 = keff*hbar/me;
mfp0 = vth0.*taoband0;
Ptrf0 = exp(-a./mfp0);
%modified band mobility at source end
taoeff0 = taoband0.*(1+a./mfp0);
miueff0 = taoeff0*e/me;
%electric field at source end
Fl0 = Id(iVd)./(W*Pmtr0.*Ptrf0.*miueff0.*ntot0.*e);
%% mobility and electric field across the channel
%define channel profiles
Vy = zeros(nL,1);
ntrap = zeros(length(Vy),1);
nband = zeros(length(Vy),1);
ncorr = zeros(length(Vy),1);
ntot = zeros(length(Vy),1);
ntot(1) = ntot0;
Ef = zeros(length(Vy),1);
tao_tc = zeros(length(Vy),1);
taoband = zeros(length(Vy),1);
miuband = zeros(length(Vy),1);
Pmtr = zeros(length(Vy),1);
mfp = zeros(length(Vy),1);
Ptrf = zeros(length(Vy),1);
taoeff = zeros(length(Vy),1);
miueff = zeros(length(Vy),1);
Fl = zeros(length(Vy),1);
Fl(1) = Fl0;
vdi = zeros(length(Vy),1);
Idf = Id(iVd);
scr = 0.05;
%for each postion in the channel (starting from source end)
for iVy = 2:length(Vy)
itnFl = 1;
crtFl = 0;
crtFl(itnFl) = 1;
Fl(iVy) = Fl(iVy-1);
Vy(iVy) = trapz(Lpo(1:iVy),Fl(1:iVy))+Vy(1);
vdi(iVy) = vdi(iVy-1);
Flagsat = 0;
while abs(crtFl(itnFl)) >= 1E-2
if Vy(iVy) <= Vov(iVd) && Flagsat == 0
Flagsat = 0;
scr = 0.1;
Fl(iVy) = Fl(iVy)+scr*(Id(iVd)-Idf)/Id(iVd)*Fl(iVy);
Vy(iVy) = trapz(Lpo(1:iVy),Fl(1:iVy));
dFldy = (Fl(iVy)-Fl(iVy-1))/dL;
ncorr(iVy) = d*Es*E0/e*dFldy;
ntot(iVy) = Cox*(Vov(iVd)-Vy(iVy))/e+ncorr(iVy);
if ntot(iVy) <= 0
Flagsat = 1;
end
end
if Vy(iVy) > Vov(iVd) || Flagsat == 1
Flagsat = 1;
scr = 0.1;
if ntot(iVy) <= 0
ntot(iVy) = ntot(iVy-1);
ncorr(iVy) = ntot(iVy)-Cox*(Vov(iVd)-Vy(iVy))/e;
elseif ntot(iVy) > 0
ntotf = Id(iVd)/W/e/vdi(iVy);
ntot(iVy) = ntot(iVy)+scr*(ntotf-ntot(iVy));
ncorr(iVy) = ntot(iVy)-Cox*(Vov(iVd)-Vy(iVy))/e;
end
Fl(iVy) = (ntot(iVy)-Cox*(Vov(iVd)-Vy(iVy))/e)/(d*Es*E0/e/dL)+Fl(iVy-1);
Vy(iVy) = trapz(Lpo(1:iVy),Fl(1:iVy));
end
%solve fermi level and Pmtr
Teff = sqrt(T^2+(gamma*e*Fl(iVy)*a/(kb)).^2);
syms En Efs
g_band = gs*gv*me/(2*pi*hbar^2);
DOS_band = g_band.*heaviside(En-Ec);
Tta = Nt/(g_band)/kb;
g_trap = Nt/(kb*Tta)*exp(En/(kb*Tta));
DOS_trap = g_trap.*heaviside(Ec-En);
DOS_tot = DOS_band + DOS_trap;
FE = 1./(1+exp((En-Efs)/(kb*Teff)));
ntotint = FE*DOS_tot;
eqn = ntot(iVy) == vpaintegral(ntotint,En,-2*e,2*e);
Efsol = double(vpasolve(eqn,Efs));
if isempty(Efsol) == 1
Ef(iVy) = 1*e;
elseif isempty(Efsol) == 0
Ef(iVy) = Efsol;
end
Eint = e*linspace(-5,5,1000);
FEint = 1./(1+exp((Eint-Ef(iVy))/(kb*Teff)));
ntrapint = FEint.*Nt/(kb*Tta).*exp(Eint/(kb*Tta)).*heaviside(Ec-Eint);
nbandint = FEint.*gs*gv*me/(2*pi*hbar^2).*heaviside(Eint-Ec);
ntrap(iVy) = trapz(Eint,ntrapint);
nband(iVy) = ntot(iVy)-ntrap(iVy);
Pmtr(iVy) = 1-ntrap(iVy)./ntot(iVy);
%% band mobility
% trapped carrier scattering
kav = sqrt(2*kb*T*me)/hbar; %non-degenerate
kf = sqrt(2*pi*nband(iVy)); %degenerate
theta = linspace(0,pi)’;
S1 = e^2*nband(iVy)/(2*kb*T*E0*Es); %non-degenerate
S2 = 2*e^2*me/(4*pi*E0*Es*hbar^2); %degenerate
S0 = S1.*heaviside(S2-S1)+S2.*heaviside(S1-S2);
keff = kav.*heaviside(S2-S1)+kf.*heaviside(S1-S2);
Ibbint = (sin(theta)).^2./(sin(theta)+S0/keff).^2;
Ibb = trapz(theta,Ibbint);
tao_tc(iVy) = 1./(e^4*(ntrap(iVy))/d*me/(8*pi*(E0*Es)^2*hbar^3*(keff).^3)*Ibb);
% PO_phonon scattering
Ediff = 1/(1/Esinf-1/Es);
Np = 1/(exp(hbar*w0/(kb*T))-1);
tao_po = 4*pi*E0*Ediff*hbar^2/(e^2*w0*me*d)/Np;
% band mobility
taoband(iVy) = 1/(1/tao_po+1/tao_tc(iVy));
miuband(iVy) = taoband(iVy)*e/me;
% calculate Ptrf
vth = hbar*keff/me;
mfp(iVy) = vth.*taoband(iVy);
Ptrf(iVy) = exp(-a./mfp(iVy));
taoeff(iVy) = taoband(iVy).*(1+a./mfp(iVy));
miueff(iVy) = taoeff(iVy)*e/me;
% carrier ensemble velocity and drain current
vsat0 = sqrt(8*hbar*w0/(3*pi*me));
vdi(iVy) = Pmtr(iVy).*Ptrf(iVy).*miueff(iVy).*Fl(iVy)./sqrt(1+(Fl(iVy).*miueff(iVy)./vsat0).^2);
Idf = W*ntot(iVy).*e.*vdi(iVy);
itnFl = itnFl+1;
crtFl(itnFl) = max(abs((Idf-Id(iVd))./Id(iVd)));
end
end
Vd(iVd) = Vy(end);
end
%% result
ntot(1) = ntot0;
nband(1) = nband0;
ntrap(1) = ntrap0;
Pmtr(1) = Pmtr0;
Ptrf(1) = Ptrf0;
miueff(1) = miueff0;
%end