%Fig4_3.m %Essential Electron Transport for Device Physics %Input: Potential profile, doping, temperature %Output: Current density as function of applied voltage bias %This program uses % - CalculateChemicalPotential.m to calculate the chemical potential Mu % for given temperature and doping concentration % - PotentialProfileWithPoissonB.m to calculate the form of the potential % profile in the charge accumulation/depletion region by solving Poisson % equation in depletion approximation % - option input parameters from file 'TunnelCurrentParameterInputB.txt' % - option input potential profile from file % 'TunnelCurrentPotentialProfileInputB.txt' clear all; clf; FS = 12; %label fontsize 18 FSN = 12; %number fontsize 16 LW = 1; %linewidth % Change default axes fonts. set(0,'DefaultAxesFontName', 'Times'); set(0,'DefaultAxesFontSize', FSN); % Change default text fonts. set(0,'DefaultTextFontname', 'Times'); set(0,'DefaultTextFontSize', FSN); %t = cputime; %start timer eye = complex(0.,1.); %square root of -1 hbar = 1.05457159e-34; %Planck's constant (Js) echarge = 1.6021764e-19; %electron charge (C) kB = 1.3806505e-23/echarge;%Boltzmann constant (eV/K) m0 = 9.109382e-31; %bare electron mass (kg) eps0 = 8.854188e-12; %permittivity of free space (F/m) VbiasMin = (pi*10^-9); %min. applied voltage bias (eV) VbiasMax = 0.40; %max. applied voltage bias (eV) VbiasPoints = 26; %number of discrete voltage bias points, 26 Vbias = linspace(VbiasMin,VbiasMax,VbiasPoints);%Voltage bias array % input parameters from file 'ParameterInput.txt': %[ParameterList]=textread('TunnelCurrentParameterInputB.txt','%*s %f %*s'); % m = ParameterList(1) * m0; %effective electron mass (kg) % eps1 = ParameterList(2) * eps0; %permittivity of GaAs 13.1 (F/m) % eps2 = ParameterList(3) * eps0; %permittivity of AlAs 10.1 (F/m) % monolayer = ParameterList(4); %thickness of one monolayer GaAs (m) % Temperature = ParameterList(5); %temperature (K) % NDopingcm = ParameterList(6) % meff 0.07 (m/m0) % eps1 13.1 (eps(GaAs)/eps0) % eps2 13.1 (eps(AlAs)/eps0) % monolayer 0.282665e-9 (m) % Temperature 300.0 (K) % NDoping 1.0 (10^18/cm^(-3)) xi = 0.36; %Al fraction in GaAs V0 = 0.84*xi; %heterostructure conduction band offset (eV) m = 0.07 * m0; %effective electron mass (kg) 0.07 eps1 = 13.1 * eps0; %permittivity of GaAs 13.1 (F/m) eps2 = 13.1 * eps0; %permittivity of AlAs 10.1 (F/m) monolayer = 0.282665e-9; %thickness of one monolayer GaAs (m) 0.282665nm Temperature = 300.0; %temperature (K) 300 NDopingcm = 1.0; %doping x10^18 cm^-3 NDoping = NDopingcm*10^24; %doping concentration in contacts (m^-3) %end reading in parameters if Temperature < 1; Temperature = 1; end %never allow T=0 K beta=1/(kB*Temperature); %inverse thermal energy (eV^-1) EFermi = ((hbar^2)*((3.0*(pi^2)*NDoping)^(2.0/3.0)))/(2.0*m*echarge); Emax = EFermi + (15.0/beta); %cutoff for energy integral (eV) nEnergy=2*VbiasPoints; %number of energy steps in integration Emin = (3*pi*10^-8); %avoid divide by zero dEnergy = Emax/nEnergy; const0 = (m*(echarge^3)*dEnergy*10^-4)/(beta*2*(pi^2)*(hbar^3));%10^-4 to convert from A/m^2 to A/cm^2 const2 = 2.0*echarge*m/(hbar^2); %chemical potential for given temperature and doping Mu = CalculateChemicalPotential(Temperature,NDoping,m,EFermi); %set up energy array and precalculate bias independent terms energyindex=[1:nEnergy]; Energy=Emin+dEnergy*double(energyindex); expa=(beta*(Mu-Energy)); const1a = 1.0+(exp(expa)); nEnergyx2=2*nEnergy; eindex=2*energyindex; %setup potential profile: %read data L(#monolayer), V(eV) from file 'PotentialProfileInput.txt' %L(#monolayer) denotes the width of barrier region in the potential %profile with a certain value of potential V(eV) %e.g.: for a single 4 monolayer thick 1 eV energy rectangular barrier, %we have 1 barrier region with L(1) = 4 monolayer and V(1) = 1 eV %[L, V] = textread('TunnelCurrentPotentialProfileInputB.txt','%f %f','headerlines',2); %L = [0.01 0.01 0.01 0.01 5 5 5 4 4 5 5 5 0.01 0.01 0.01 0.01 0.01]; %V = [0 0 0 0 V0 V0 0 0 0 0 V0 V0 0 0 0 0 0]; L = [0.01 5 5 5 4 4 5 5 5 0.01]; V = [0 0 0 V0 V0 V0 V0 0 0 0]; L = monolayer*L;%thickness of barrier in the potential profile (m) NBarrier = length(L);%NBarrier = # of barriers in the potential profile Lx(1) = 0.0; for i = 2 : 1 : NBarrier+1 Lx(i) = Lx(i-1) + L(i-1); end d2 = Lx(NBarrier+1); %total thickness of potential profile (m) d1Max = -(eps1*d2)/(2.0*eps2)+sqrt(VbiasMax/echarge*eps1/NDoping+... (eps1*d2)/(2.0*eps2)*(eps1*d2)/(2.0*eps2)); %Width of charge accumulation/depletion region %width of depletion region for maximum value of voltage bias, %depending on NDoping xmin = -(d1Max+2.0e-9); %lower cutoff in x-direction, depending on NDoping xmax = d2+d1Max+2.0e-9; %upper cutoff in x-direction, depending on NDoping xpoints = round((xmax-xmin)*1.0e10)*2; %# of points in x-direction for propagation matrix method dx = (xmax-xmin)/xpoints; E2=Energy.'*ones(1,xpoints); for i = 1 : 1 : xpoints %set up x-axis x(i) = double(i-1) * dx + xmin; end % fprintf(1,'\n'); % fprintf(1,'Voltage bias (V) Current (A/cm^2)\n'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %outer loop calculates potential profile at each voltage bias Vbias(j) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Current=zeros(1,VbiasPoints);%initialize t = cputime; %start timer for j = 1:1:VbiasPoints %calculate potential for each Vbias(j) Vpot = PotentialProfileWithPoissonB(xpoints, x, dx, Lx, NBarrier, V, Vbias(j), NDoping, eps1, eps2); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %loop to integrate over energy starting from Emin and stepping by dEnergy %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% expb=(beta*(Mu-Energy-Vbias(j))); const1b =1.0+(exp(expb)); const1ab=const1a./const1b; const1=log(const1ab); V2=ones(nEnergy,1)*Vpot; bigP2=speye(nEnergyx2); k2=sqrt(const2*(E2-V2)); for n=1:1:xpoints-1 k2tmp1=exp(-eye*k2(:,n)*dx); k2tmp2=exp(eye*k2(:,n)*dx); p11=(1+k2(:,n+1)./k2(:,n)).*k2tmp1; p12=(1-k2(:,n+1)./k2(:,n)).*k2tmp1; p21=(1-k2(:,n+1)./k2(:,n)).*k2tmp2; p22=(1+k2(:,n+1)./k2(:,n)).*k2tmp2; p2=sparse(eindex-1,eindex-1,p11,nEnergyx2,nEnergyx2); p2=p2+sparse(eindex-1,eindex,p12,nEnergyx2,nEnergyx2); p2=p2+sparse(eindex,eindex-1,p21,nEnergyx2,nEnergyx2); p2=p2+sparse(eindex,eindex,p22,nEnergyx2,nEnergyx2); bigP2=0.5*bigP2*p2; end Trans=diag(bigP2); Trans=Trans(eindex-1); Trans=abs(1./Trans).^2; Trans=Trans.*const1'; Current(j)=const0*sum(Trans); % fprintf(1,' %1.3f %5.9f \n', Vbias(j), Current(j)); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %end of outer Vbias loop %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% run_time = cputime-t %CPU run time %plot results CurrentMax=max(Current); ttl=(['\rmFig4.3, \itn\rm=',num2str(NDopingcm,'%5.1f'),... '\times10^{18} cm^{-3}, \itT\rm=',... num2str(Temperature,'%5.0f'),' K, \mu=',num2str(Mu*1000,'%5.1f'),... ' meV, \itm\rm^*_e=',num2str(m/m0,'%7.2f'),'\timesm_0, \itV\rm_0=',... num2str(V0,'%7.2f'),' eV']); figure(1); plot(Vbias,Current,'b','LineWidth',LW); xlabel('Applied voltage, \itV\rm_{bias} (V)','FontSize',FS); ylabel('Current density, \itJ\rm_{tot} (A/cm^{2})','FontSize',FS) axis([0, VbiasMax, 0, 1.1*CurrentMax]); grid on title(ttl,'FontSize',FS); figure(2); plot(Vbias,log10(Current),'b','LineWidth',LW); xlabel('Applied voltage, \itV\rm_{bias} (V)','FontSize',FS); ylabel('log_{10} (\itJ\rm_{tot}) (A cm^{-2})','FontSize',FS); axis([0, VbiasMax, log10(Current(2)), log10(1.1*CurrentMax)]); grid on title(ttl,'FontSize',FS); figure(3); plot(x*1e9,Vpot,'b','LineWidth',LW); xlabel('Position, \itx\rm (nm)','FontSize',FS); ylabel('Potential, \itV\rm (eV)','FontSize',FS); axis([min(x*1e9),max(x*1e9), 1.1*min(Vpot),1.1*max(Vpot)]); title(ttl,'FontSize',FS);