Fall 2017: ECE 5390 / MSE 5472: Quantum Transport in Electron Devices and Novel Materials

The Ballistic Transistor

h=((6.63*(10^(-34)))/(2*π)) ;  (* Reduced Planck's Constant *)
q=1.6*(10^(-19));      (* Electron charge, Coulomb *)
m0=9.1*(10^(-31));   (* Free electron mass, kg *)
mt=0.2*m0;         (* Transverse effective mass, Silicon *)
mGaN=0.2*m0;     (* Electron effective mass of GaN *)
k=1.38*(10^(-23));    (* Boltzmann Constant, J/K *)
eps0=8.85*(10^(-12));   (* Permitivitty of vacuum, F/m *)
epsOX=3.9*eps0;   (* Dielectric constant of SiO2 *)
epsAlOX=9.0*eps0;   (* Dielectric constant of SiO2 *)
epsAlN=9.0*eps0;  (* Dielectric constant of AlN *)
Mv=2.5;  (* Valley degeneracy of Silicon *)
Cb[epsb_,tb_]:=((epsb*eps0)/(tb*(10^(-9))))  (* in F/(m^2), epsb is relative dielectric constant, tb barrier thickness in nm *)
Cq[meff_,gv_]:= (((q^2)*2*gv*meff*m0)/(2*π*(h^2))) ; (* Quantum capacitance, in F/(m^2), gv is valley degeneracy *)
Vth[T_]:=((k*T)/q); (* thermal voltage in Volts *)
vd[VD_,T_]:=(VD/(Vth[T]));  (* dimensionless drain voltage *)
J02d[meff_,gv_,T_]:=((q*gv*(Sqrt[2*meff*m0])*(((k*T))^(3/2)))/((π^2)*(h^2))); (* spin degeneracy is gs=2, current in A/m *)
nsvg[tb_,epsb_,meff_,gv_,T_,Vgs_]:=FindRoot[Exp[((q*(ns*(10^16)))/(Vth[T]))*((1/(Cq[meff,gv]))+(1/(Cb[epsb,tb])))]-Exp[((q*(ns*(10^16)))/(Vth[T]))*((1/(Cb[epsb,tb])))]==Exp[(Vgs/(Vth[T]))],{ns,15}][[1]][[2]]
Quiet[LogPlot[{nsvg[1,4,0.2,1,300,Vgs],((10^(-16))/(q*((1/(Cq[0.2,1]))+(1/(Cb[4,1])))))*Vgs,(10^(-16))*((Cq[0.2,1]*Vth[300])/q)*Exp[(Vgs/(Vth[300]))]},{Vgs,-0.8,1},Frame->True,PlotStyle->{{Black,Thick},{Black,Dashed},{Black,Dashed}},FrameLabel->{"(Subscript[V, gs])-(Subscript[V, T])(Volt)","(Subscript[n, s])((10^12)/(cm^2))"},PlotRange->{(10^(-8)),(10^2)},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]]

Gamma[(1/2)+1]
(∫_(-(π/2))^(π/2))Cos[θ]\[DifferentialD]θ
(∫_0^(π/2))Sin[θ]Cos[θ]\[DifferentialD]θ

(* Fermi Dirac Integral of order j=1/2 *)
F[u_]:=(2/(Sqrt[π]))*NIntegrate[((Sqrt[y])/(1+Exp[y-u])),{y,0,100}] (*Fermi Dirac Integral of order j=1/2 *)
Quiet[LogPlot[{F[u],-PolyLog[(3/2),-(E^u)],Exp[u],((u^(3/2))/(3/2))},{u,-10,10},PlotRange->{(10^(-6)),(10^2)},Frame->True,PlotStyle->{{Black,Thick},{Red,Dashed},{Black,Dashed},{Black,Dashed}},FrameLabel->{"u","Fermi Dirac integrals"},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic] ](* checking values of the Fermi-Dirac integral *)

Quiet[Plot[{F[u],-PolyLog[(3/2),-(E^u)],Exp[u],((u^(3/2))/(3/2))},{u,-10,10},PlotRange->{(10^(-6)),(10^2)},Frame->True,PlotStyle->{{Black,Thick},{Red,Dashed},{Black,Dashed},{Black,Dashed}},FrameLabel->{"u","Fermi Dirac integrals"},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic] ](* checking values of the Fermi-Dirac integral *)

(*Fermi Dirac Integral of order j=1/2 *)
ηs[Vgs_,VT_,tb_,epsb_,meff_,gv_,T_,VD_]:=Log[Exp[((q*(nsvg[tb,epsb,meff,gv,T,Vgs]*(10^16)))/(Cq[meff,gv]*Vth[T]))]-1];
(* ηg[VG_,VT_,tb_,epsb_,meff_,gv_,T_]:=((Cb[epsb,tb])/(Cb[epsb,tb]+Cq[meff,gv]))*((VG-VT)/(Vth[T]));  (* dimensionless gate voltage parameter *)
ηs[VG_,VT_,tb_,epsb_,meff_,gv_,T_,VD_]:=Log[√((((1+Exp[vd[VD,T]]))^2)+4*Exp[vd[VD,T]]*(Exp[2*ηg[VG,VT,tb,epsb,meff,gv,T]]-1))-(1+Exp[vd[VD,T]])]-Log[2.0];  (* dimensionless source voltage parameter *) *)
J2dballistic[VG_,VT_,tb_,epsb_,meff_,gv_,T_,VD_]:=(10^(-3))*J02d[meff,gv,T]*(F[ηs[VG,VT,tb,epsb,meff,gv,T,VD]]-F[ηs[VG,VT,tb,epsb,meff,gv,T,VD]-vd[VD,T]]);  (* Ballistic current per unit width in mA/micron *)
Quiet[LogPlot[{J2dballistic[Vg,0,2,10,0.2,2.5,300,0.05],J2dballistic[Vg,0,2,10,0.2,2.5,300,0.5]},{Vg,-1,1},Frame->True,PlotRange->{(10^(-10)),10},FrameLabel->{"(Subscript[V, gs])-(Subscript[V, T])","(Subscript[I, d])/W, mA/micron"},PlotLabel->"Ballistic Si MOSFET, T=300K",PlotStyle->{Thick,Dashed},BaseStyle->{FontSize->15},AspectRatio->GoldenRatio,GridLines->Automatic]]
(* Ballistic Silicon MOSFET at T=300K *)
(* J2dballistic[VG,VT,tb,epsb,meff,gv,T,VD] *)
Quiet[Plot[{J2dballistic[1.0,0,2,10,0.2,2.5,300,VD],J2dballistic[0.75,0,2,10,0.2,2.5,300,VD],J2dballistic[0.5,0,2,10,0.2,2.5,300,VD],J2dballistic[0.25,0,2,10,0.2,2.5,300,VD],J2dballistic[0.001,0,2,10,0.2,2.5,300,VD]},{VD,0,0.5},Frame->True,PlotRange->{-0.1,5},FrameLabel->{"(Subscript[V, ds])","(Subscript[I, d])/W, mA/micron"},PlotLabel->"Ballistic Si MOSFET, T=300K",PlotStyle->Thick,BaseStyle->{FontSize->15}] ](* valley degeneracy gv=2.5, meff=0.2*m0, gate sweeps Vgs-Vt=1.25 V to 0.0 V in steps of 0.25 V, tb=2 nm, dielectric epsb=10 *)

Fermi’s Golden Rule

2-level system excited by a step perturbation

ConSq[W12_,δe_]:=((W12^2)/((((((δe)/2)))^2)+(W12^2)));
Ω[W12_,δe_]:=(Sqrt[(((((δe)/2)))^2)+(W12^2)]);
c1sq[t_,W12_,δe_]:=1-ConSq[W12,δe]*((Sin[Ω[W12,δe]*t])^2);
c2sq[t_,W12_,δe_]:=ConSq[W12,δe]*((Sin[Ω[W12,δe]*t])^2);
Manipulate[Plot[{c1sq[t,W12,δe],c2sq[t,W12,δe]},{t,0,100},Frame->True,PlotRange->{0,1},PlotStyle->Thick,FrameLabel->{"time","Probability"},BaseStyle->{FontSize->20}],{W12,0.02,1},{δe,0,1}]
Manipulate[Plot[ConSq[W12,δe],{δe,-1,1},Frame->True,PlotRange->{0,1},PlotStyle->Thick,FrameLabel->{"δE","(C^2)"},BaseStyle->{FontSize->20}],{W12,0.02,1}]
Manipulate[Plot[ConSq[W12,δe],{W12,0,5},Frame->True,PlotRange->{0,1},PlotStyle->Thick,FrameLabel->{"(Subscript[W, 12])","(C^2)"},BaseStyle->{FontSize->20}],{δe,0.02,1}]

Fermi’s Golden Rule

Manipulate[Plot[(((((Sin[((ω*t)/2)])/(ω/2))))^2),{ω,-10,10},Frame->True,FrameLabel->{"ω=(((Subscript[ϵ, n]) - (Subscript[ϵ, 0]))/ℏ)","(((Subscript[Pr, |0⟩ -> |n⟩])[t])/(|⟨n | V | 0⟩((|)^2)))"},PlotStyle->Thick,BaseStyle->{FontSize->15},PlotRange->{-0.5,20},Filling->Bottom],{t,0,10*((2π)/5)}]

Manipulate[Plot[(((((Sin[((ω*t)/2)])/(ω/2))))^2),{ω,-10,10},Frame->True,FrameLabel->{"ω=(((Subscript[ϵ, n]) - (Subscript[ϵ, 0]))/ℏ)","(((Subscript[Pr, |0⟩ -> |n⟩])[t])/(|⟨n | V | 0⟩((|)^2)))"},PlotStyle->Thick,BaseStyle->{FontSize->15},PlotRange->{-0.01,1},Filling->Bottom],{t,0,1*((2π)/5)}]

(∫_(-∞)^∞)(((((Sin[((ω*t)/(2ℏ))])/(ω/2))))^2)\[DifferentialD]ω

fg[x_,η_]:=((2*η)/((x^2)+(η^2)));
Manipulate[Plot[fg[x,η],{x,-10,10},PlotRange->{0,10},Frame->True,PlotStyle->Thick,Filling->Axis,BaseStyle->{FontSize->15},FrameLabel->{"x","((2*η)/((x^2) + (η^2)))"}],{η,10,0}]

Manipulate[Plot[fg[x,η],{η,-1,1},PlotRange->{-50,50},Frame->True,PlotStyle->Thick,Filling->Axis,BaseStyle->{FontSize->15},FrameLabel->{"η","((2*η)/((x^2) + (η^2)))"}],{x,0.1,0.000000001}]

TransRate[t_,eta_,won_,w_]:=Exp[2*hbeta*t]*(((2*eta)/((((won+w))^2)+(((eta))^2)))+((2*eta)/((((won-w))^2)+(((eta))^2))))*(1-Cos[2*w*t])+ 2*Sin[2*w*t]*(((won+eta)/((((won+w))^2)+(((eta))^2)))+((won-w)/((((won-w))^2)+(((eta))^2))))

(∫_(-∞)^t)Exp[(I*(Subscript[ω, nm])+η)*tp]\[DifferentialD]tp//FullSimplify

(∫_(-∞)^t)(∫_(-∞)^tp)Exp[(I*(Subscript[ω, nm])+η)*tp]*Exp[(I*(Subscript[ω, m0])+η)*tpp]\[DifferentialD]tpp\[DifferentialD]tp//FullSimplify

(∫_(-∞)^t)(∫_(-∞)^tp)(∫_(-∞)^tpp)Exp[(I*(Subscript[ω, nm])+η)*tp]*Exp[(I*(Subscript[ω, ml])+η)*tpp]*Exp[(I*(Subscript[ω, l0])+η)*tppp]\[DifferentialD]tppp\[DifferentialD]tpp\[DifferentialD]tp//FullSimplify

Boltzmann Transport Equation

(* The De-Broglie wavelength *)
λdB[meff_,T_]:=((6.63*(10^(-34)))/(Sqrt[2*meff*9.1*(10^(-31))*0.026*((T/300))*1.6*(10^(-19))]))*(10^9);
Plot[{λdB[0.05,T],λdB[0.2,T],λdB[0.4,T]},{T,100,400},FrameLabel->{"Temperature (K)","de Broglie wavelength (nm)"},GridLines->Automatic,Frame->True,PlotStyle->Thick,BaseStyle->{FontSize->15},PlotRange->{0,50}]

(* The Fermi-Dirac Distribution Function *)
fd0[kx_,ky_,Ec_,Ef_,T_]:=(1/(1+Exp[((Ec+((((1.055*(10^(-34))))^2)/(2*0.2*9.1*(10^(-31))*1.6*(10^(-19))))((((kx*(10^8)))^2)+(((ky*(10^8)))^2))-Ef)/(0.026*((T/300))))]));  (* 2D Fermi-Dirac distribution *)
Manipulate[Plot[fd0[kx,0,0,0.1,T],{kx,-10,10},PlotRange->{0,1.2},Frame->True,PlotStyle->Thick,FrameLabel->{"(Subscript[k, x])((10^8)/m)","f((Subscript[k, x]),(Subscript[k, y])=0) "},BaseStyle->{FontSize->15}],{T,500,10}]
Manipulate[Plot3D[fd0[kx,ky,0,0.1,T],{kx,-20,20},{ky,-20,20},PlotRange->{0,1.2}],{T,500,10}]

(* Perturbed Occupation Function *)
f0[kx_,ky_,Ec_,Ef_,T_]:=(1/(1+Exp[((Ec+((((1.055*(10^(-34))))^2)/(2*0.2*9.1*(10^(-31))*1.6*(10^(-19))))(((((kx-2)*(10^8)))^2)+(((ky*(10^8)))^2))-Ef)/(0.026*((T/300))))]));  (* perturbed 2D distribution function *)
Manipulate[Plot[f0[kx,0,0,0.1,T],{kx,-10,10},PlotRange->{0,1.2},Frame->True,PlotStyle->Thick,FrameLabel->{"(Subscript[k, x])((10^8)/m)","f((Subscript[k, x]),(Subscript[k, y])=0) "},BaseStyle->{FontSize->15}],{T,500,10}]
Manipulate[Plot3D[f0[kx,ky,0,0.1,T],{kx,-20,20},{ky,-20,20},PlotRange->{0,1.2}],{T,500,10}]

(* Transient Fermi occupation function *)
ftrans[kx_,ky_,Ec_,Ef_,T_,t_,taum_]:=fd0[kx,ky,Ec,Ef,T]+(f0[kx,ky,Ec,Ef,T]-fd0[kx,ky,Ec,Ef,T])*Exp[-(t/taum)];
Manipulate[Plot[ftrans[kx,0,0,0.1,T,t,0.2],{kx,-10,10},PlotRange->{0,1.2},Frame->True,PlotStyle->Thick,FrameLabel->{"(Subscript[k, x])((10^8)/m)","f((Subscript[k, x]),(Subscript[k, y])=0) "},BaseStyle->{FontSize->15}],{t,0,1},{T,500,10}]
Manipulate[Plot3D[ftrans[kx,ky,0,0.1,T,t,0.2],{kx,-20,20},{ky,-20,20},PlotRange->{0,1.2}],{t,0,1},{T,500,10}]

(* Derivative of the Fermi function showing Dirac-delta behavior *)
fd0[kx_,ky_,Ec_,Ef_,T_]:=(1/(1+Exp[((Ec+((((1.055*(10^(-34))))^2)/(2*0.2*9.1*(10^(-31))*1.6*(10^(-19))))((((kx*(10^8)))^2)+(((ky*(10^8)))^2))-Ef)/(0.026*((T/300))))]));  (* 2D Fermi-Dirac distribution *)
df0dkx[kx_,ky_,Ec_,Ef_,T_]=D[fd0[kx,ky,Ec,Ef,T],kx];
Manipulate[Plot[df0dkx[kx,0,0,0.1,T],{kx,-10,10},PlotRange->{-2,2},Frame->True,PlotStyle->Thick,FrameLabel->{"(Subscript[k, x])((10^8)/m)","df/dkx((Subscript[k, x]),(Subscript[k, y])=0) "},BaseStyle->{FontSize->15}],{T,500,10}]

Mobility and Scattering

Neutral impurity

(* Seeger Eq 6.1.9 *)
muNeutral[Nx_,meff_,epsDc_]:=((1.44*(10^22))/(Nx*(10^17)))*(meff/epsDc); (* Neutral impurity scattering limited mobility in (cm^2)/Vs, Nx is in (10^17)/(cm^3) *)
LogLogPlot[muNeutral[Nx,0.2,9],{Nx,0.1,10},Frame->True,AspectRatio->1,PlotStyle->Thick,FrameLabel->{"Nx ((10^17)/(cm^3))","Mobility ((cm^2)/Vs)"},BaseStyle->{FontSize->15}]

Charged (Ionized) Impurity

(* Seeger Eq 6.3.23 *)
(∫_0^∞)Exp[-(r/LD)]*Sin[q*r]\[DifferentialD]r
βBH[epsDc_,T_,meff_,n_]:=(Sqrt[(epsDc/16)])*((T/100))*(Sqrt[meff])*(Sqrt[((2.08*(10^18))/(n*(10^17)))]);
muIonizedImp[Nd_,epsDc_,T_,meff_,n_]:=((3.68*(10^20))/(Nd*(10^17)))*((((epsDc/16)))^2)*((((T/100)))^(3/2))*(1/(((Sqrt[meff])/(Log[10.0]))*(Log[1+((βBH[epsDc,T,meff,n])^2)]-(((βBH[epsDc,T,meff,n])^2)/(1+((βBH[epsDc,T,meff,n])^2)))))); (* Ionized impurity scattering limited mobility in (cm^2)/Vs, Nd is the ionized impurity density and n is the free carrier concentration, both in (10^17)/(cm^3) *)
LogLogPlot[{muIonizedImp[Nd,16,300,1,Nd],muIonizedImp[Nd,16,77,1,Nd]},{Nd,(10^(-4)),(10^1)},Frame->True,PlotRange->{(10^3),(10^8)},AspectRatio->1,PlotStyle->Thick,FrameLabel->{"(Subscript[N, D]) ((10^17)/(cm^3))","Mobility ((cm^2)/Vs)"},BaseStyle->{FontSize->15}] (* Fig 6.4 Seeger *)
LogLogPlot[{muIonizedImp[Nd,9,300,0.2,Nd],muIonizedImp[Nd,9,77,0.2,Nd]},{Nd,(10^(-4)),(10^1)},PlotRange->{(10^3),(10^8)},Frame->True,AspectRatio->1,PlotStyle->Thick,FrameLabel->{"(Subscript[N, D]) ((10^17)/(cm^3))","Mobility ((cm^2)/Vs)"},BaseStyle->{FontSize->15}] (* For GaN *)

Acoustic Phonon (Deformation potential)

Phon[Num_,x_,t_]:=(∑_(n=1)^Num)n*Re[Exp[I*(x-n*t)]]
Manipulate[Plot[Phon[Num,x,t],{x,-5,5},PlotStyle->Thick,Frame->True,PlotRange->{-20,20}],{t,0,100},{Num,1,5,1}]

(* Seeger Eq. 6.4.18 *)
muacdp[T_,ac_,meff_,cl_]:=(10^4)*((2(Sqrt[2*π]))/3)*(((1.6*(10^(-19)))*(((1.0552*(10^(-34))))^4)*(cl*(10^12)*((10^(-5))/(10^(-4)))))/((((meff*9.1*(10^(-31))))^(5/2))*(((ac*1.6*(10^(-19))))^2)*(((1.6*(10^(-19))*0.026*(T/300)))^(3/2)))); (* Acoustic phonon scattering limited mobility in (cm^2)/Vs, ac is deformation potential in eV, cl is in (10^12)dynes/(cm^2), 1 dyne = (10^(-5)) Newton, *)
(* muacdp[100,9.5,0.2,1.56] (* calibration for Ge, get same result as Seeger *) *)
LogLogPlot[muacdp[T,8.4,0.2,2.65],{T,10,300},Frame->True,AspectRatio->1,PlotStyle->Thick,FrameLabel->{"T (K)","Mobility ((cm^2)/Vs)"},BaseStyle->{FontSize->15}] (* For GaN *)

Acoustic Phonon (Piezoelectric)

K[epsDc_,epz_,cl_]:=(((((epz*1.6*(10^(-19))*(10^14)))^2)/(cl*(10^12)*((10^(-5))/(10^(-4)))))/((epsDc*8.5*(10^(-12)))+((((epz*1.6*(10^(-19))*(10^14)))^2)/(cl*(10^12)*((10^(-5))/(10^(-4))))))) ;(* epz is the peizoelectric coefficient in units of (10^12) electron charges/(cm^2) *)
K[10,1,1.8]

Optical Phonon (Deformation Potential)

fint[z_]:=NIntegrate[(((y^(3/2))*Exp[-2*z*y])/((Sqrt[y+1])+Exp[2*z]*Re[(Sqrt[y-1])])),{y,0,1000}];
Plot[fint[z],{z,0,1}]

Optical Phonon (Polar)

(* Seeger Eq 6.12.13 *)
α[Eop_,meff_,epsInf_,epsDc_]:=397.4*(Sqrt[(meff/(((Eop/26)*300)))])*((1/epsInf)-(1/epsDc));  (* Dimensionless Frohlich coupling constant, optical phonon energy expressed in meV *)
θD[Eop_]:=((Eop/26)*300); (* Debye temperature *)
mupop3d[T_,Eop_,meff_,epsInf_,epsDc_]:=(10^4)*((1.6*(10^(-19)))/(2*(9.1*(10^(-31))*meff)*α[Eop,meff,epsInf,epsDc]*((((10^(-3))*1.6*(10^(-19))*Eop)/(1.0552*(10^(-34)))))))*Exp[((θD[Eop])/T)];  (* Polar optical phonon scattering limited mobility in (cm^2)/Vs *)
LogLogPlot[{mupop3d[T,36,0.067,10.9,12.53] (*GaAs*),mupop3d[T,92,0.2,5,9] (*GaN*), mupop3d[T,55,2,50,300](*STO*)},{T,1,300},PlotRange->{10,(10^7)},Frame->True,AspectRatio->1,PlotStyle->Thick,FrameLabel->{"T (K)","Mobility ((cm^2)/Vs)"},BaseStyle->{FontSize->15}]

Dislocation

(* Seeger Eq 6.15.6 *)
LD[epsDc_,T_,Nd_]:=(Sqrt[((epsDc*8.85*(10^(-12))*1.6*(10^(-19))*26*(10^(-3))*(T/300))/((((1.6*(10^(-19))))^2)*Nd*(10^22)))]);  (* Non degenerate Debye length in m *)
mudisl3d[Ndis_,a_,T_,meff_,Nd_,epsDc_,f_]:=(10^3)*((30*(Sqrt[2*π])*(((epsDc*8.85*(10^(-12))))^2)*(((a*(10^(-10))))^2))/((((1.6*(10^(-19))))^3)*(f^2)*LD[epsDc,T,Nd]*(Sqrt[meff*9.1*(10^(-31))])))*((((1.6*(10^(-19))*26*(10^(-3))*(T/300)))^(3/2))/(Ndis*(10^12))); (*Dislocation scattering limited mobility in (cm^2)/Vs, Ndis is in (10^8)/(cm^2), a is in Angstroms, Nd is in (10^16)/(cm^3) *)
LogLogPlot[{mudisl3d[1,5,T,0.2,30,10,1],mudisl3d[10,5,T,0.2,30,10,1],mudisl3d[100,5,T,0.2,30,10,1]},{T,1,300},Frame->True,AspectRatio->1,PlotStyle->Thick,FrameLabel->{"T (K)","Mobility ((cm^2)/Vs)"},BaseStyle->{FontSize->15}]

Alloy Disorder

All combined

Manipulate[LogLogPlot[{mupop3d[T,92,0.2,5,9],mudisl3d[Ndis,5,T,0.2,30,10,1],muacdp[T,8.4,0.2,2.65],muIonizedImp[Nd,9,300,0.2,Nd],muNeutral[Nx,0.2,9],1/((1/(mupop3d[T,92,0.2,5,9]))+(1/(mudisl3d[Ndis,5,T,0.2,30,10,1]))+(1/(muacdp[T,8.4,0.2,2.65]))+(1/(muIonizedImp[Nd,9,300,0.2,Nd]))+(1/(muNeutral[Nx,0.2,9])))},{T,1,300},Frame->True,PlotRange->{(10^1),(10^7)},AspectRatio->1,PlotStyle->Thick,FrameLabel->{"T (K)","Mobility ((cm^2)/Vs)"},BaseStyle->{FontSize->15}],{Ndis,0.01,100},{Nd,0.1,50},{Nx,0.01,5}] (* Plotting for GaN, be aware of the approximations such as non-degenerate carriers, etc *)

The Long-Channel Transistor

(* Define constants *)
h=((6.63*(10^(-34)))/(2*π)) ;  (* Planck's Constant *)
q=1.6*(10^(-19));      (* Electron charge, Coulomb *)
m0=9.1*(10^(-31));   (* Free electron mass, kg *)
mt=0.2*m0;         (* Transverse effective mass, Silicon *)
mGaN=0.2*m0;     (* Electron effective mass of GaN *)
kb=1.38*(10^(-23));    (* Boltzmann Constant, J/K *)
eps0=8.85*(10^(-12));   (* Permitivitty of vacuum, F/m *)
epsOX=3.9*eps0;   (* Dielectric constant of SiO2 *)
epsAlOX=9.0*eps0;   (* Dielectric constant of SiO2 *)
epsAlN=9.0*eps0;  (* Dielectric constant of AlN *)
Mv=2.5;  (* Valley degeneracy of Silicon *)

ψB[T_,ni_,Na_]:=((kb*T)/q)*Log[((Na*(10^23))/(ni*(10^16)))];
ψsRoots[Vgs_,Vfb_,epsS_,epsB_,tox_,T_,ni_,Na_,V_]:=FindRoot[Vgs-Vfb-ψs==((Sqrt[2*epsS*eps0*kb*T*(Na*(10^23))])/((eps0*epsB)/(tox*(10^(-9)))))*(Sqrt[((q*ψs)/(kb*T))+(((((ni*(10^16))/(Na*(10^23)))))^2)*Exp[((q*(ψs-V))/(kb*T))]]),{ψs,1}][[1]][[2]]

Quiet[Plot[{2*ψB[300,1,1]+V,ψsRoots[1,0,9,3.9,10,300,1,1,V],ψsRoots[2,0,9,3.9,10,300,1,1,V],ψsRoots[3,0,9,3.9,10,300,1,1,V]},{V,0,2.5},Frame->True,PlotStyle->Thick,FrameLabel->{"Electron quasi Fermi level V (Volt)","Surface Potential (Subscript[ψ, s])(Volt)"},BaseStyle->{FontSize->15},PlotRange->{0,3}]]

Manipulate[Quiet[Plot[{2*ψB[300,1,1]+V,ψsRoots[Vgs,0,9,3.9,10,300,1,1,V]},{V,0,2.5},Frame->True,PlotStyle->Thick,FrameLabel->{"Electron quasi Fermi level V (Volt)","Surface Potential (Subscript[ψ, s])(Volt)"},BaseStyle->{FontSize->15},PlotRange->{0,3}]],{Vgs,0,4}]

IdsByW[μ_,L_,Vgs_,Vfb_,epsS_,epsB_,tox_,T_,ni_,Na_,V_]:=(μ*(10^(-4)))*(1/((L*(10^(-6)))))*(((eps0*epsB)/(tox*(10^(-9))))*(Vgs-Vfb+((kb*T)/q))*ψsRoots[Vgs,Vfb,epsS,epsB,tox,T,ni,Na,V]-(1/2)*((eps0*epsB)/(tox*(10^(-9))))*((ψsRoots[Vgs,Vfb,epsS,epsB,tox,T,ni,Na,V])^2)-(2/3)*(Sqrt[2*eps0*epsS*q*(Na*(10^23))])*((ψsRoots[Vgs,Vfb,epsS,epsB,tox,T,ni,Na,V])^(3/2))+((kb*T)/q)*(Sqrt[2*eps0*epsS*q*(Na*(10^23))*ψsRoots[Vgs,Vfb,epsS,epsB,tox,T,ni,Na,V]]));
IdsTot[μ_,L_,Vgs_,Vfb_,epsS_,epsB_,tox_,T_,ni_,Na_,Vds_]:=IdsByW[μ,L,Vgs,Vfb,epsS,epsB,tox,T,ni,Na,Vds]-IdsByW[μ,L,Vgs,Vfb,epsS,epsB,tox,T,ni,Na,0];
Quiet[Plot[{IdsTot[100,1,0,0,9,3.9,10,300,1,1,Vds],IdsTot[100,1,1,0,9,3.9,10,300,1,1,Vds],IdsTot[100,1,2,0,9,3.9,10,300,1,1,Vds],IdsTot[100,1,3,0,9,3.9,10,300,1,1,Vds],IdsTot[100,1,4,0,9,3.9,10,300,1,1,Vds]},{Vds,0,5},Frame->True,PlotStyle->Thick,FrameLabel->{"(Subscript[V, ds]) (Volt)","(Subscript[I, d])/W (A/m)"},BaseStyle->{FontSize->15},PlotRange->{0,100}]]

Quiet[LogPlot[{IdsTot[100,1,Vgs,0,9,3.9,10,300,1,1,0.1],IdsTot[100,1,Vgs,0,9,3.9,10,300,1,1,1]},{Vgs,0,2},Frame->True,PlotStyle->Thick,FrameLabel->{"(Subscript[V, gs]) (Volt)","Log[(Subscript[I, d])/W]"},BaseStyle->{FontSize->15},AspectRatio->1,PlotRange->{(10^(-10)),10}]]