%Fig3_7.m %Essential Electron Transport for Device Physics %resonant tunneling through 3 barriers 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); hbar_txt=['\fontname{MT Extra}h\fontname{Arial}']; t=cputime; nbarrier=3; %number of potential barriers bx=0.4e-9; %barrier width (m) wx=0.6e-9; %well width (m) V0=0.5; %barrier energy (eV) N=140+1; %number of samples of potential for j=1:1:20 %set up position and potential array dL(j)=wx/5; V(j)=0; dL(j+20)=bx/20; V(j+20)=V0; %barrier 1 dL(j+40)=wx/20; V(j+40)=0; dL(j+60)=bx/20; V(j+60)=V0; %barrier 2 dL(j+80)=wx/20; V(j+80)=0; dL(j+100)=bx/20; V(j+100)=V0; %barrier 3 dL(j+120)=wx/5; V(j+120)=0; end dL(N)=wx/20; V(N)=0; xp(1)=dL(1); for j=2:N;xp(j)=xp(j-1)+dL(j);end; %set up x-plot axis Emin=pi*1e-5; %add (pi*1.0e-5) to energy avoid divide by zero Emax=0.8; %maximum particle energy (eV) npoints=1500; %number of points in energy plot dE=Emax/npoints; %energy increment (eV) hbar=1.0545715e-34; %Planck's constant (J s) eye=complex(0.,1.); %square root of -1 m0=9.109382e-31; %bare electron mass (kg) meff=1.0; %effective electron mass / m0 m=meff*m0; %effective electron mass (kg) echarge=1.6021764e-19; %electron charge (C) % vectorized form of outer for-loop below, computing all energylevels % simultaneously energyindex=[1:npoints]; E=Emin+dE*energyindex; E2=E.'*ones(1,N); V2=ones(npoints,1)*V; npointsx2=2*npoints; bigP2=speye(npointsx2); k2=sqrt(2*echarge*m*(E2-V2))/hbar; for n=1:(N-1) p11=0.5*(1+k2(:,n+1)./k2(:,n)).*exp(-eye*k2(:,n)*dL(n)); p12=0.5*(1-k2(:,n+1)./k2(:,n)).*exp(-eye*k2(:,n)*dL(n)); p21=0.5*(1-k2(:,n+1)./k2(:,n)).*exp(eye*k2(:,n)*dL(n)); p22=0.5*(1+k2(:,n+1)./k2(:,n)).*exp(eye*k2(:,n)*dL(n)); p2=sparse(2*energyindex-1,2*energyindex-1,p11,npointsx2,npointsx2); p2=p2+sparse(2*energyindex-1,2*energyindex,p12,npointsx2,npointsx2); p2=p2+sparse(2*energyindex,2*energyindex-1,p21,npointsx2,npointsx2); p2=p2+sparse(2*energyindex,2*energyindex,p22,npointsx2,npointsx2); bigP2=bigP2*p2; end Trans2=diag(bigP2); Trans2=Trans2(2*energyindex-1); Trans2=abs(1./Trans2).^2; %************************************************************************** %Plot potential and transmission coefficient %************************************************************************** figure(1); %plot potential subplot(1,2,1) plot(xp*1e9,V,'LineWidth',LW) axis([0,xp(end)*1e9,0,Emax]) xlabel('Position, x (nm)') ylabel('Potential energy, V (eV)') ttl=['\rmFig3.7, m^*_e=',... num2str(meff),'\timesm_0, L_b=',num2str(bx*1e9),... ' nm, L_w=',num2str(wx*1e9),' nm']; title(ttl); subplot(1,2,2),plot(Trans2,E,'LineWidth',LW); axis([0,1,0,Emax]); xlabel('Transmission coefficient'); ylabel('Energy, E (eV)'); figure(2); %plot transmission coefficient subplot(1,2,1),plot(xp*1e9,V,'LineWidth',LW); axis([0,xp(end)*1e9,0,Emax]); xlabel('Position, x (nm)'); ylabel('Potential energy, V (eV)'); title(ttl); subplot(1,2,2);plot(-log(Trans2),E,'LineWidth',LW);axis([0,20,0,Emax]); xlabel('-ln trans. coeff.');ylabel('Energy, E (eV)'); %************************************************************************** %Calculate electron wave function %************************************************************************** E1=0.3218; %fixed electron energy (eV), 0.3218 0.4016 bigP=[1,0;0,1]; %default value of matrix bigP for i=1:N k(i)=sqrt(2*echarge*m*(E1-V(i)))/hbar;%wave vector at energy E1 end % vectorized construction of k k2=sqrt(2*echarge*m*(E1-V))/hbar; for n=1:(N-1) %multiply out propagation matrix p(1,1)=0.5*(1+k(n+1)/k(n))*exp(-eye*k(n)*dL(n)); p(1,2)=0.5*(1-k(n+1)/k(n))*exp(-eye*k(n)*dL(n)); p(2,1)=0.5*(1-k(n+1)/k(n))*exp(eye*k(n)*dL(n)); p(2,2)=0.5*(1+k(n+1)/k(n))*exp(eye*k(n)*dL(n)); bigP=bigP*p; end C(1)=1/bigP(1,1); %transmission amplitude B(1)=bigP(2,1)*C(1); %reflection amplitude A(1)=1; %incident amplitude bigP=[1,0;0,1]; %default value of matrix bigP for n=1:(N-1) %multiply out propagation matrix p(1,1)=0.5*(1+k(n+1)/k(n))*exp(-eye*k(n)*dL(n)); p(1,2)=0.5*(1-k(n+1)/k(n))*exp(-eye*k(n)*dL(n)); p(2,1)=0.5*(1-k(n+1)/k(n))*exp(eye*k(n)*dL(n)); p(2,2)=0.5*(1+k(n+1)/k(n))*exp(eye*k(n)*dL(n)); bigP=bigP*p; Prop=inv(bigP)*[A(1);B(1)]; C(n)=Prop(1); D(n)=Prop(2); A(n+1)=C(n); B(n+1)=D(n); end for n=1:(N-1) psi(n)=(C(n)+D(n)); %calculate wave function end psi(N)=psi(N-1); %dummy value for N-th value of wave function psi2=conj(psi).*psi; %calculate wave function squared %************************************************************************** %Plot modulus of wave function squared %************************************************************************** figure(3); plot(xp*1e9,psi2,'LineWidth',LW); xlabel('Position, x (nm)'),ylabel('|\psi(x)|^2'); ttl2=['\rm',num2str(nbarrier),'-barrier, m^*_e=',... num2str(meff),' \times m_0, L_b=',num2str(bx*1e9),... ' nm, L_w=',num2str(wx*1e9),' nm, E=',num2str(E1),' eV']; title(ttl2); %************************************************************************** %Plot real and imaginary part of wave function %************************************************************************** figure(4); plot(xp*1e9,real(psi),'b','LineWidth',LW); hold on plot(xp*1e9,imag(psi),'r','LineWidth',LW); xlabel('Position, x (nm)'),ylabel('Real (blue), Imag (red), (\psi(x))'); title(ttl2); hold off e=cputime-t; e