%Fig3_13b.m %Essential Electron Transport for Device Physics %Calculate inelastic transmission through an arbitrary %potential containing an Einstein phonon at one location. The temperature %is 0K, so that the incident electron cannot absorb a phonon. clear all clf FS = 12; %label FontSize 18 FSN = 12; %number FontSize 16 LW = 1; %LineWidth % Change default axes fonts. set(0,'DefaultAxesFontName', 'Times New Roman'); set(0,'DefaultAxesFontSize', FSN); % Change default text fonts. set(0,'DefaultTextFontname', 'Times New Roman'); set(0,'DefaultTextFontSize', FSN); hbar_txt=['\fontname{MT Extra}h\fontname{Times New Roman}']; %******************* % Input Parameters * %******************* % Phonon Parameters %omega=0.01, g1=0.07, meff=0.07, Lb=2, LW=2.95+0.01+2.95, V0=0.3, Emax=30 %omega=0.01, g1=0.07, meff=0.07, Lb=1.5, LW=2.95+0.01+2.95, V0=0.3, Emax=30 %%omega=0.01, g1=0.22, meff=1.0, Lb=0.4, LW=0.295+0.01+0.295, V0=1, Emax=150 omega = 0.036; % Phonon Energy (eV) 0.036 0.01 g0 = 0.0*1e-9; % Static Delta Barrier Energy, units of eVnm 0.2 0.0 0.0 g1 = 0.080*1e-9; % units of eVm 0.12 0.08 0.008 %g1c is 0.14 eV nm in GaAs % Electron-Phonon Coupling Constant (eV) 0.07 0.15 0.22 NumPhonons = 8; % Number of Phonons Included in Calculation 5 % Energy Parameters Eplot = 1.3774; % Energy Value for Plotting of Wavefunction EMin = pi*1e-5; % Minimum Injection Energy (eV) 0.999*omega+pi*1e-5 EMax = 3.5*omega;% Maximum Injection Energy (eV)3.5 1.001*omega NumE = 500; % Number of Samples in Energy Range 1000 Vbias = 0; % Voltage Bias (eV) % Spatial Parameters Numx = 100; % Number x-axis points for potential and wave function plots 500 % Material Parameters % Barrier Parameters Potential = [ 0.0 0.0 0.0 0.0 0.0 0.0 0.0]; % Potential Energy of Each Segment (eV) Length = 1e-9*[ 1.5 1.5 1.5 1.5 1.5 1.5 1.5]; % Length of Each Segment (m) Phonon = [ 0.0 0.0 0.0 1.0 0.0 0.0 0.0]; % 1 Indicates Which Barrier Phonon is Located at %************************* % Constants and Matrices * %************************* % Physical Constants hbar = 1.05457148e-34; % Planck's constant (Js) echarge = 1.60217646e-19; % Electron charge (C) m0 = 9.11e-31; % Electron mass (kg) meff = 0.07; % Effective electron mass mass=m0*meff; %convert g to units of eV %gconv=sqrt(meff*m0*echarge/(2*hbar)); % gconv=sqrt(2*meff*m0*echarge)/hbar; % g0=g0*gconv; % g1=g1*gconv; % Include Zero Phonon NumPhonons = NumPhonons + 1; % Initialize Arrays Transmission = zeros(NumE,1); % Transmission Array ElasticTrans = zeros(NumE,1); % Elastic Only Transmission Array TransCoeffs = zeros(NumE,NumPhonons + 1); % Transmission Coefficients Array %****************** % Build Potential * %****************** x = linspace(0,sum(Length),Numx); dx = x(2) - x(1); xPhonon = zeros(1,Numx); V = zeros(1,Numx); BiasSlope = Vbias/(sum(Length) - Length(1) - Length(end)); for j = 1:Numx for k = 1:length(Length) if (x(j) >= sum(Length(1:k-1))) && (x(j) < sum(Length(1:k))) V(j) = Potential(k); if (x(j) > Length(1)) && (x(j) <= (sum(Length) - Length(end))) V(j) = V(j) - (x(j) - Length(1))*BiasSlope; elseif x(j) > (sum(Length) - Length(end)) V(j) = V(j) - Vbias; end end end end V(end) = Potential(end) - Vbias; % Find phonon location PhononLocation = 0; for j = 1:length(Length) if Phonon(j) == 1 break else PhononLocation = PhononLocation + Length(j); end end [z,PhononIndex] = min(abs(x - PhononLocation)); xPhonon(PhononIndex) = 1; clear z PhononIndex PhononLocation %****************************************************** % Calculate Transmission and Reflection Due to System * %****************************************************** % No Phonons if g1=0 if g1 == 0 NumPhonons = 1; end %Create Injection Energy Array E = linspace(EMin,EMax,NumE); % Find E value closest to Eplot Ediff = (E - Eplot).^2; [z,EplotIndex] = min(Ediff); clear z % Loop Over Injection Energies for j = 1:NumE P = eye(2*NumPhonons); % Reset P ElasticP = eye(2); for k = 1:Numx-1 % Loop through x-axis % Create k arrays for n = 1:NumPhonons+2 k1(n) = sqrt(2*mass*echarge*(E(j)-(n-2)*omega - V(k)))/hbar; k2(n) = sqrt(2*mass*echarge*(E(j)-(n-2)*omega - V(k+1)))/hbar; end if xPhonon(k) == 1 % Create Inelastic Step Matrix Pstep = zeros(2*NumPhonons,2*NumPhonons+4); for n = 1:NumPhonons AlphaN = 1i*mass*echarge*g1*sqrt(n-1)*ones(1,2)/(k1(n+1)*hbar^2); BetaN = 1i*mass*echarge*g1*sqrt(n)*ones(1,2)/(k1(n+1)*hbar^2); Pstep(2*n-1,2*n-1:2*n) = AlphaN; Pstep(2*n,2*n-1:2*n) = -AlphaN; Pstep(2*n-1,2*n+3:2*n+4) = BetaN; Pstep(2*n,2*n+3:2*n+4) = -BetaN; Pstep(2*n-1,2*n+1) = (k1(n+1) + k2(n+1))/(2*k1(n+1)) + (1i*mass*echarge*g0)/(hbar^2*k1(n+1)); Pstep(2*n-1,2*n+2) = (k1(n+1) - k2(n+1))/(2*k1(n+1)) + (1i*mass*echarge*g0)/(hbar^2*k1(n+1)); Pstep(2*n,2*n+1) = (k1(n+1) - k2(n+1))/(2*k1(n+1)) - (1i*mass*echarge*g0)/(hbar^2*k1(n+1)); Pstep(2*n,2*n+2) = (k1(n+1) + k2(n+1))/(2*k1(n+1)) - (1i*mass*echarge*g0)/(hbar^2*k1(n+1)); end Pstep = Pstep(:,3:2*NumPhonons+2); else % Create Elastic Step Matrix Pstep = zeros(2*NumPhonons); for n = 1:NumPhonons Pstep(2*n-1,2*n-1) = (k1(n+1) + k2(n+1))/(2*k1(n+1)); Pstep(2*n,2*n) = Pstep(2*n-1,2*n-1); Pstep(2*n,2*n-1) = (k1(n+1) - k2(n+1))/(2*k1(n+1)); Pstep(2*n-1,2*n) = Pstep(2*n,2*n-1); end end ElasticStep = zeros(2); % Calculate Elastic Step Matrix if xPhonon(k) == 0 ElasticStep(1,1) = k1(2) + k2(2); ElasticStep(2,2) = ElasticStep(1,1); ElasticStep(1,2) = k1(2) - k2(2); ElasticStep(2,1) = ElasticStep(1,2); ElasticStep = ElasticStep/(2*k1(2)); else ElasticStep(1,1) = k1(2) + k2(2) + 1i*2*echarge*mass*g0/(hbar^2); ElasticStep(2,2) = ElasticStep(1,1); ElasticStep(1,2) = k1(2) - k2(2) + 1i*2*echarge*mass*g0/(hbar^2); ElasticStep(2,1) = ElasticStep(1,2); ElasticStep = ElasticStep/(2*k1(2)); end % Create Propagation Matrix Pprop = zeros(2*NumPhonons); for n = 1:NumPhonons Pprop(2*n-1,2*n-1) = exp(-1i*k2(n+1)*dx); Pprop(2*n,2*n) = exp(1i*k2(n+1)*dx); end % Create Inelastic Propagation Matrix ElasticProp = zeros(2); ElasticProp(1,1) = exp(-1i*k2(2)*dx); ElasticProp(2,2) = exp(1i*k2(2)*dx); % Multiply Out Transfer Matrices P = P*Pstep*Pprop; ElasticP = ElasticP*ElasticStep*ElasticProp; end % Move bn's Into T Matrix for n = 1:NumPhonons P(:,2*n) = zeros(2*NumPhonons,1); P(2*n,2*n) = -1; end % Create A Matrix A = zeros(2*NumPhonons,1); A(1,1) = 1; % Solve System T = P\A;%note that P\A is faster that inv(P)*A % Store Largest n That has a Propagating Mode Leaving the Barrier nMax = 0; for n = 1:NumPhonons if E(j) > (n-1)*omega + V(end) nMax = n; end end % Calculate k values for transmission kin = sqrt(2*mass*echarge*(E(j) - V(1)))/hbar; kout = zeros(1,nMax); for n = 1:nMax kout(n) = sqrt(2*mass*echarge*(E(j) - (n-1)*omega - V(end)))/hbar; end % Add Transmission Coefficients for Propogating Modes for n = 1:nMax Transmission(j) = Transmission(j) + abs(kout(n)*T(2*n-1)*conj(T(2*n-1))/kin); TransCoeffs(j,n+1) = abs(kout(n)*T(2*n-1)*conj(T(2*n-1))/kin); end TransCoeffs(j,1) = Transmission(j); % Save data for wavefunction plots if j == EplotIndex TPlot = T; end % Solve elastic system ElasticTrans(j) = 1/abs(ElasticP(1,1))^2; end %************************** % Calculate Wavefunctions * %************************** % Eliminate b values from T matrix for n = 1:NumPhonons TPlot(2*n) = 0; end % Initialize arrays and matrices for wavefunction calculation Psi = zeros(2*NumPhonons,Numx); Psi(:,end) = TPlot; PsiPlot = zeros(NumPhonons,Numx); % Back Calculate Wavefunction Values for j = 1:Numx - 1 k = Numx - j; % Create k arrays for n = 1:NumPhonons+2 k1(n) = sqrt(2*mass*echarge*(E(EplotIndex)-(n-2)*omega - V(k)))/hbar; k2(n) = sqrt(2*mass*echarge*(E(EplotIndex)-(n-2)*omega - V(k+1)))/hbar; end if xPhonon(k) == 1 % Create Inelastic Step Matrix Pstep = zeros(2*NumPhonons,2*NumPhonons+4); for n = 1:NumPhonons AlphaN = 1i*echarge*g1*sqrt(n-1)*ones(1,2)/(k1(n+1)); BetaN = 1i*echarge*g1*sqrt(n)*ones(1,2)/(k1(n+1)); Pstep(2*n-1,2*n-1:2*n) = AlphaN; Pstep(2*n,2*n-1:2*n) = -AlphaN; Pstep(2*n-1,2*n+3:2*n+4) = BetaN; Pstep(2*n,2*n+3:2*n+4) = -BetaN; Pstep(2*n-1,2*n+1) = ((k1(n+1) + k2(n+1))/(2*k1(n+1))) + (1i*mass*echarge*g0)/(hbar^2*k1(n+1)); Pstep(2*n-1,2*n+2) = ((k1(n+1) - k2(n+1))/(2*k1(n+1))) + (1i*mass*echarge*g0)/(hbar^2*k1(n+1)); Pstep(2*n,2*n+1) = ((k1(n+1) - k2(n+1))/(2*k1(n+1))) - (1i*mass*echarge*g0)/(hbar^2*k1(n+1)); Pstep(2*n,2*n+2) = ((k1(n+1) + k2(n+1))/(2*k1(n+1))) - (1i*mass*echarge*g0)/(hbar^2*k1(n+1)); end Pstep = Pstep(:,3:2*NumPhonons+2); else % Create Elastic Step Matrix Pstep = zeros(2*NumPhonons); for n = 1:NumPhonons Pstep(2*n-1,2*n-1) = (k1(n+1) + k2(n+1))/(2*k1(n+1)); Pstep(2*n,2*n) = Pstep(2*n-1,2*n-1); Pstep(2*n,2*n-1) = (k1(n+1) - k2(n+1))/(2*k1(n+1)); Pstep(2*n-1,2*n) = Pstep(2*n,2*n-1); end end % Create Propagation Matrix Pprop = zeros(2*NumPhonons); for n = 1:NumPhonons Pprop(2*n-1,2*n-1) = exp(-1i*k2(n+1)*dx); Pprop(2*n,2*n) = exp(1i*k2(n+1)*dx); end Psi(:,k) = Pstep*Pprop*Psi(:,k+1); end % Build and Normalize Wavefunctions for k = 1:Numx for n = 1:NumPhonons PsiPlot(n,k) = Psi(2*n-1,k) + Psi(2*n,k); end end % Normalize Wavefunctions and Find Zero Values for Plots WaveMax = zeros(1,NumPhonons); for n = 1:NumPhonons if max(PsiPlot(n,:)) ~= 0 WaveMax(n) = max(PsiPlot(n,:)); PsiPlot(n,:) = 0.9*PsiPlot(n,:)*omega/(2*WaveMax(n)); end PsiZero(n) = E(EplotIndex) - (n-1)*omega; end PsiMagPlot = 3.5*PsiPlot.*conj(PsiPlot); %*************** % Plot Results * %*************** % Create logarithmic transmissions TransLog = -log(Transmission); % Calculate Transmission Derivatives TransDiff = diff(Transmission); Trans2Diff = diff(TransDiff); % Set delta magnitude if max(V) == zeros(1,Numx) DeltaMax = 1; else DeltaMax = 1.2*max(V)*sign(max(V)); end ttl=(['\rmFig3.13b', hbar_txt,'\omega=',... num2str(omega),' eV, \itm\rm_{eff}=',num2str(meff),... ', \itg\rm_0=',num2str(g0*1e9),' eV nm, \itg\rm_1=',... num2str(g1*1e9),' eV nm, \itN\rm_p=',num2str(NumPhonons-1)]); %************************************************************************** % Plot transmission by itself figure(1) clf %subplot(1,2,1) plot(E,ElasticTrans,'k--','Linewidth',LW); hold on plot(E,Transmission,'Linewidth',LW) axis([EMin EMax min(Transmission) 1.002]) %axis([EMin EMax 0 1.1]) title(ttl) xlabel(['Energy, \itE\rm (eV)']) ylabel('Transmission') legend('Elastic','Inelastic','Location','Southeast') grid on %************************************************************************** % Plot transmission and potential for comparison figure(2) clf subplot(1,3,1) plot(Transmission,E/omega,'Linewidth',LW) title(ttl) xlabel('Transmission') ylabel(['Normalized energy, \itE\rm/ ',hbar_txt,'\omega']) axis([0 1 1.2*min(EMin,min(V))/omega max(EMax/omega,max(xPhonon))]) grid on subplot(1,3,2) plot(x*1e9,DeltaMax*xPhonon/omega,'r','Linewidth',LW) hold on plot(x*1e9,V/omega,'k','Linewidth',LW) hold off %title('Potential (delta (red))') xlabel('\itx\rm, (nm)') ylabel(['\itV\rm/ ',hbar_txt,'\omega']) axis([x(1)*1e9 x(end)*1e9 1.2*min(EMin,min(V))/omega max(EMax/omega,max(xPhonon))]) grid on subplot(1,3,3) plot(TransLog,E/omega,'Linewidth',LW) %title('ln T') xlabel('-ln Transmission') ylabel(['Normalized energy, \itE\rm/ ',hbar_txt,'\omega']) axis([0 max(TransLog) 1.2*min(EMin,min(V))/omega max(EMax/omega,max(xPhonon))]) grid on %************************************************************************** % Plot Channel Coefficients figure(3) clf plot(E,TransCoeffs,'Linewidth',LW) title(ttl) xlabel(['Energy, \itE\rm (eV)']) ylabel('Transmission coefficients') axis([EMin EMax 0 1]) grid on % Make Normalization Text if NumPhonons == 1 Normalizations = []; else Normalizations = ['Multipliers for wavefunctions are ']; Normalizations2= ['Multipliers for densities are ']; for n = 1:NumPhonons Normalizations =[Normalizations num2str(abs(WaveMax(1)/WaveMax(n))) ', ']; Normalizations2=[Normalizations2 num2str(abs(WaveMax(1)/WaveMax(n))^2) ', ']; end end %************************************************************************** return %************************************************************************** % Plot Transmission Derivatives figure(4) clf subplot(1,2,1) plot(E(1:end-1)/omega,TransDiff,'Linewidth',2) title('First Derivative of Transmission') xlabel('Normalized Energy (E/\Omega)') ylabel('dT/dE') xlim([E(1)/omega E(end-1)/omega]) grid on subplot(1,2,2) plot(E(1:end-2)/omega,Trans2Diff,'Linewidth',2) title('Second Derivative of Transmission') xlabel('Normalized Energy (E/\Omega)') ylabel('d^2T/dE^2') xlim([E(1)/omega E(end-1)/omega]) grid on % Plot Wavefunctions figure(5) clf subplot(1,3,1) hold on for n = 1:NumPhonons plot(x*1e9,real(PsiPlot(n,:) + PsiZero(n)),'b') plot(x*1e9,imag(PsiPlot(n,:) + i*PsiZero(n)),'r') plot(x*1e9,PsiZero(n)*ones(1,Numx),'g') end hold off title({['\rmWavefunctions for inelastically scattered electron with injection energy of ',num2str(E(EplotIndex)),' eV'];['Wavefunctions offset so that the zero level is E_{inj} - n\Omega and normalized to prevent overlap'];Normalizations}) xlabel('\itx\rm (nm)') ylabel('Wavefunction, \psi (m^{-1/2})') axis([x(1)*1e9 x(end)*1e9 min((E(EplotIndex)-(NumPhonons-0.5)*omega),0) 1.2*max((E(EplotIndex)+omega/2),max(V))]) subplot(1,3,2) plot(x*1e9,DeltaMax*xPhonon,'r') hold on plot(x*1e9,V,'k') hold off title(['\rmBackground Potential (black) and delta phonon (red) with a coupling constant of ',num2str(g1),' eV']) xlabel('\itx\rm (nm)') ylabel('Potential Energy (eV)') axis([x(1)*1e9 x(end)*1e9 min((E(EplotIndex)-(NumPhonons-0.5)*omega),0) 1.2*max((E(EplotIndex)+omega/2),max(V))]) subplot(1,3,3) hold on for n = 1:NumPhonons plot(x*1e9,real(PsiMagPlot(n,:)) + PsiZero(n)) plot(x*1e9,PsiZero(n)*ones(1,Numx),'g') end hold off title({['Probability density for inelastically scattered electron with injection energy of ',num2str(E(EplotIndex)),' eV'];['Wavefunctions offset so that the zero level is E_{inj} - n\Omega and normalized to prevent overlap'];Normalizations2}) xlabel('\itx\rm (nm)') ylabel('Probability Density, |\psi|^2 (m^{-1})') axis([x(1)*1e9 x(end)*1e9 min((E(EplotIndex)-(NumPhonons-0.5)*omega),0) 1.2*max((E(EplotIndex)+omega/2),max(V))])