%Fig3_19.m %Essential Electron Transport for Device Physics clear all clc clf warning off 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 = .036; % Phonon energy (eV) g0 = 0*1e-9; % Static delta barrier energy (eVm) g1 = 0.1*1e-9; % Electron-phonon coupling constant (eVm) 0.1 0.5 NumPhonons = 3; % Number of phonons included in calculation NumPhononsPlot = 4; % +1 Number of channels shown in movie when inelastic scattering occurs (must be >= 1) % Bias Parameters Bias = 0.0; % Potential bias applied to device (V) VShift = 0.0; % Offset for entire potential profile (V) % Wavenumber Parameters kCenter = 2e8; % Center wavenumber of gaussian (m^-1) (2.1432e8 = 25 meV) sigmak = 0.5e8; % Broadening Full-Width-Half-Maximum (eV) (sigmax = 1/sigmak) 0.5 NumGauss = 200; % Number of Gaussian frequency components % Spatial Parameters Numx = 170; % Number of points in x-axis for potential and wave function plots 151 % Material Parameters meff = .07; % Effective mass coefficient NDoping = 1e18; % Carrier concentration (cm^-3) epsr1 = 13.1; % Relative permittivity of outer regions epsr2 = 13.1; % Relative permittivity of barrier regions % Temporal Parameters tMin = 0; % Start time RunTimeMultiplier = 2; % Run time of movie in terms of the time it takes for the center of the wave packet to traverse the screen % e.g. Setting equal to 1 will end movie with center of wave packet at right edge of screen Numt = 350; % Number of samples of time range 250 % Barrier Parameters Potential = [0.0 -0.075 -0.075 -0.025]; % Potential in each section Length = [70 15 15 70]*1e-9; % Length of each section Phonon = [0 0 1 0]; % Phonon locations. A '1' puts a delta barrier at the left edge of the section. Put a '2' in last column to put delta at right edge of domain. %************************* % Constants and matrices * %************************* % Universal Constants hbar = 1.0545715e-34; % Reduced Planck's Constant (J s) m0 = 9.109382e-31; % Bare Electron Mass (kg) echarge = 1.6021764e-19; % Electron Charge (C) eps0 = 8.854188e-12; % Permittivity of Free Space (F/m) % Modify Parameters m = meff*m0; % Effective electron mass (kg) eps1 = epsr1*eps0; % Permittivity of Ohmic Contact eps2 = epsr2*eps0; % Permittivity of Ohmic Contact NDoping = NDoping*1e6; % Modify doping concentration to m^-3 NumPhonons = NumPhonons + 1; % Add zero phonon % Calculate Parameters tMax = RunTimeMultiplier*m*sum(Length)/(hbar*kCenter); ECenter = hbar^2*kCenter^2/(2*m*echarge) ESigmak = hbar^2*sigmak^2/(2*m*echarge) % If any parameter indicates no phonons, clear all phonon related parameters if sum(Phonon) == 0 || g1 == 0 || omega == 0 || NumPhonons == 1 g1 = 0; NumPhonons = 1; omega = 0; end % Initialize Matrices kAmps = zeros(1, NumGauss); psi = zeros(2*NumPhonons, Numx, NumGauss); Psi = zeros(Numt, Numx, NumPhonons); %********************************************* % Set up arrays and matrices for calculation * %********************************************* % Create k array k = linspace(kCenter - 3*sigmak, kCenter + 3*sigmak, NumGauss); % Warn user if negative k values present if k(1) < 0 disp('Warning: Negative values of k are in Gaussian wave packet and will be ignored') end % Create energy array Energy = hbar^2*k.*k/(2*m*echarge); % Calculate energy amplitudes kAmps = sqrt(sigmak/sqrt(pi))*exp(-((kCenter - k).*(kCenter - k))/(2*sigmak^2))/(sigmak*sqrt(2*pi)); % Create Time Array Time = linspace(tMin, tMax, Numt); % Make new barrier length array NumBarrier = length(Length); Length = [0 Length]; Steps = Length*triu(ones(NumBarrier+1)); % Calculate Poisson parameters d2 = Steps(NumBarrier+1); %Total width of potential profile (m) d1Max = -(eps1*d2)/(2.0*eps2)+sqrt(Bias*eps1/(echarge*NDoping)+(eps1*d2)^2/(2.0*eps2)^2); % Width of accumulation/depletion region % Set x-axis limits xMin = -(d1Max + 2.0e-9); %lower cutoff in x-direction, depends on NDoping xMax = d2 + d1Max + 2.0e-9;%upper cutoff in x-direction, depends on NDoping % Create x-axis x = linspace(xMin,xMax,Numx); dx = mean(diff(x)); % Locate phonon in space xPhonon = FindPhonons(Steps, Phonon, x); % Calculate potential profile V = PotentialProfileWithPoisson(x, Steps, Potential, VShift, Bias, NDoping, eps1, eps2); % Fix end of potential if no bias if Bias == 0 for xIndex = 1:Numx if V(Numx + 1 - xIndex) == Potential(end) break else V(Numx + 1 - xIndex) = Potential(end); end end end %************************** % Calculate Wavefunctions * %************************** for kIndex = 1:NumGauss % Calculate output wave function amplitudes [a, TransCoeffs, WaveAmps] = InelasticTransmissionSingleE(V, xPhonon, Energy(kIndex), dx, m, omega, g0, g1, NumPhonons); a = InelasticWavefunction(WaveAmps, V, xPhonon, Energy(kIndex), dx, m, omega, g0, g1, NumPhonons); % Get total wave functions [psi(:,:,kIndex)] = InelasticWavefunction(WaveAmps, V, xPhonon, Energy(kIndex), dx, m, omega, g0, g1, NumPhonons); % Clear variables for next run clear a TransCoeffs WaveAmps end %**************************** % Calculate Current Density * %**************************** for PhononIndex = 1:NumPhonons for TimeIndex = 1:Numt for kIndex = 1:NumGauss % Assemble Total Wavefunctions Psi(TimeIndex, :, PhononIndex) = Psi(TimeIndex, :, PhononIndex) + kAmps(kIndex)*(psi(2*PhononIndex - 1, :, kIndex) + psi(2*PhononIndex, :, kIndex))*exp(-i*(Energy(kIndex) - (PhononIndex - 1)*omega)*echarge*Time(TimeIndex)/hbar); end end end PsiMag = Psi.*conj(Psi); %******************************************** % Calculate Width and Amplitude Versus Time * %******************************************** sigmax = sqrt(sigmak^-2 + hbar^2*sigmak^2*Time.*Time/(4*m^2)); %*************** % Create Movie * %*************** for PhononIndex = 1:NumPhonons if omega == 0 Maxes(PhononIndex) = max(max(PsiMag(:, :, PhononIndex)))/0.5; else Maxes(PhononIndex) = max(max(PsiMag(:, :, PhononIndex)))/(0.75*omega); end end % Set up movie parameters avifps = ceil(Numt/10); if omega == 0 MovieMax = max(ECenter + 0.5,1.2*max(V)); MovieMin = min(0, 1.2*min(V)); else MovieMax = max(ECenter + omega,1.2*max(V)); MovieMin = min(1.2*min(Potential), 1.2*(ECenter - (NumPhononsPlot - 1)*omega)); end % Open movie file avimov = VideoWriter('Test1.avi'); open(avimov); %avimov = avifile('Test1.avi','fps',avifps); % Record movie %******************************************************************* ttl=(['\rmFig3.19, ',hbar_txt,'\omega=',... num2str(omega),' eV, \itm\rm^*_e=',num2str(meff),... '\times\itm\rm_0, \itg\rm_0=',num2str(g0*1e9),' eV nm, \itg\rm_1=',... num2str(g1*1e9),' eV nm, \itn\rm_{max}=',num2str(NumPhonons-1)]); figure(1) for TimeIndex = 1:Numt clf plot(x*1e9, V, 'k','LineWidth',LW) title(ttl); hold on if (g0 ~= 0 || g1 ~= 0) && sum(Phonon) ~= 0 DeltaLoc = find(xPhonon == 1); for DeltaIndex = 1:sum(xPhonon) line([x(DeltaLoc(DeltaIndex))*1e9 x(DeltaLoc(DeltaIndex))*1e9], [V(DeltaLoc(DeltaIndex)) .5], 'Color', 'r','LineWidth',2) end end if omega == 0 plot(x*1e9, ECenter + PsiMag(TimeIndex, :, PhononIndex)/Maxes(PhononIndex)) else for PhononIndex = 1:NumPhononsPlot plot(x*1e9, ECenter - (PhononIndex - 1)*omega + PsiMag(TimeIndex, :, PhononIndex)/Maxes(PhononIndex),'LineWidth',LW) end end hold off xlim([x(1)*1e9 x(end)*1e9]) ylim([MovieMin MovieMax]) axis square xlabel(['Position, \itx\rm (nm)']) ylabel('|\psi(\itx\rm)|^2, \itV\rm(\itx\rm) (eV)') Mov(TimeIndex) = getframe; writeVideo(avimov,Mov(TimeIndex)); %avimov = addframe(avimov,Mov(TimeIndex)); end % Close movie file %avimov = close(avimov); close(avimov); %********************************************************************** %snapshot at Numt/16 figure(2) plot(x*1e9, V, 'k','LineWidth',LW) title(ttl); hold on if (g0 ~= 0 || g1 ~= 0) && sum(Phonon) ~= 0 DeltaLoc = find(xPhonon == 1); for DeltaIndex = 1:sum(xPhonon) line([x(DeltaLoc(DeltaIndex))*1e9 x(DeltaLoc(DeltaIndex))*1e9], [V(DeltaLoc(DeltaIndex)) .5], 'Color', 'r','LineWidth',2) end end if omega == 0 plot(x*1e9, ECenter + PsiMag(ceil(Numt/6), :, PhononIndex)/Maxes(PhononIndex)) else for PhononIndex = 1:NumPhononsPlot plot(x*1e9, ECenter - (PhononIndex - 1)*omega + PsiMag(ceil(Numt/6), :, PhononIndex)/Maxes(PhononIndex),'LineWidth',LW) end end hold off xlim([x(1)*1e9 x(end)*1e9]) ylim([MovieMin MovieMax]) axis square xlabel(['Position, \itx\rm (nm)']) ylabel('|\psi(\itx\rm)|^2, \itV\rm(\itx\rm) (eV)') %******************************************************************** %space time plots %******************************************************************** xp=linspace(0,sum(Length)*1e9,Numx); yp=linspace(tMin,tMax*1e12,Numt); ttl=(['\rmFig3.19, ',hbar_txt,'\omega = ',... num2str(omega),' eV, \itg\rm_0 = ',num2str(g0*1e9),' eV nm, \itg\rm_1 = ',... num2str(g1*1e9),' eV nm, \itn\rm = ',num2str(0)]); figure(3) pcolor(xp,yp,(PsiMag(:,:,1))); colormap(jet);shading interp;%winter cool autumn jet gray copper bone xlabel(['Position, \itx\rm (nm)']) ylabel('Time, \itt\rm (ps)') title(ttl); ttl=(['\rmFig3.19, ',hbar_txt,'\omega = ',... num2str(omega),' eV, \itg\rm_0 = ',num2str(g0*1e9),' eV nm, \itg\rm_1 = ',... num2str(g1*1e9),' eV nm, \itn\rm = ',num2str(1)]); figure(4) pcolor(xp,yp,(PsiMag(:,:,2))); colormap(jet);shading interp; xlabel(['Position, \itx\rm (nm)']) ylabel('Time, \itt\rm (ps)') title(ttl); ttl=(['\rmFig3.19, ',hbar_txt,'\omega = ',... num2str(omega),' eV, \itg\rm_0 = ',num2str(g0*1e9),' eV nm, \itg\rm_1 = ',... num2str(g1*1e9),' eV nm, \itn\rm = ',num2str(2)]); figure(5) pcolor(xp,yp,(PsiMag(:,:,3))); colormap(jet);shading interp; xlabel(['Position, \itx\rm (nm)']) ylabel('Time, \itt\rm (ps)') title(ttl); ttl=(['\rmFig3.19, ',hbar_txt,'\omega = ',... num2str(omega),' eV, \itg\rm_0 = ',num2str(g0*1e9),' eV nm, \itg\rm_1 = ',... num2str(g1*1e9),' eV nm, \itn\rm = ',num2str(3)]); figure(6) pcolor(xp,yp,(PsiMag(:,:,4))); colormap(jet);shading interp; xlabel(['Position, \itx\rm (nm)']) ylabel('Time, \itt\rm (ps)') title(ttl); %located to the far right figure(7) plot(yp,log10(PsiMag(:,ceil(0.9*Numx),1)),'k','LineWidth',LW); hold on plot(yp,log10(PsiMag(:,ceil(0.9*Numx),2)),'b','LineWidth',LW); plot(yp,log10(PsiMag(:,ceil(0.9*Numx),3)),'g','LineWidth',LW); plot(yp,log10(PsiMag(:,ceil(0.9*Numx),4)),'r','LineWidth',LW); hold off xlabel('Time, \itt\rm (ps)') ylabel(['log_{10}(|\psi(\itx\rm, \itt\rm)|^2)']) title(ttl); figure(8) plot(yp,(PsiMag(:,ceil(0.9*Numx),1))./max(PsiMag(:,ceil(0.9*Numx),1)),'k','LineWidth',LW); hold on plot(yp,(PsiMag(:,ceil(0.9*Numx),2))./max(PsiMag(:,ceil(0.9*Numx),2)),'b','LineWidth',LW); plot(yp,(PsiMag(:,ceil(0.9*Numx),3))./max(PsiMag(:,ceil(0.9*Numx),3)),'g','LineWidth',LW); plot(yp,(PsiMag(:,ceil(0.9*Numx),4))./max(PsiMag(:,ceil(0.9*Numx),4)),'r','LineWidth',LW); hold off xlabel('Time, \itt\rm (ps)') ylabel(['Normalized |\psi(\itx\rm, \itt\rm)|^2']) title(ttl); %located at phonon figure(9) plot(yp,log10(PsiMag(:,ceil(0.5*Numx),1)),'k','LineWidth',LW); hold on plot(yp,log10(PsiMag(:,ceil(0.5*Numx),2)),'b','LineWidth',LW); plot(yp,log10(PsiMag(:,ceil(0.5*Numx),3)),'g','LineWidth',LW); plot(yp,log10(PsiMag(:,ceil(0.5*Numx),4)),'r','LineWidth',LW); hold off xlabel('Time, \itt\rm (ps)') ylabel(['log_{10}(|\psi(\itx\rm, \itt\rm)|^2)']) title(ttl); figure(10) plot(yp,(PsiMag(:,ceil(0.5*Numx),1))./max(PsiMag(:,ceil(0.5*Numx),1)),'k','LineWidth',LW); hold on plot(yp,(PsiMag(:,ceil(0.5*Numx),2))./max(PsiMag(:,ceil(0.5*Numx),2)),'b','LineWidth',LW); plot(yp,(PsiMag(:,ceil(0.5*Numx),3))./max(PsiMag(:,ceil(0.5*Numx),3)),'g','LineWidth',LW); plot(yp,(PsiMag(:,ceil(0.5*Numx),4))./max(PsiMag(:,ceil(0.5*Numx),4)),'r','LineWidth',LW); hold off xlabel('Time, \itt\rm (ps)') ylabel(['Normalized |\psi(\itx\rm, \itt\rm)|^2']) axis([0 max(yp) 0 1.1]); title(ttl);