%Fig7_14a.m %Essential Electron Transport for Device Physics %Inelastic electron scattering rate at finite temperature 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); colormap(jet); hbar = 1.05457159e-34; % Planck's constant (Js) nDoping1 = 1e17; % n-type doping concentration (cm^-3) 1e18 nDoping = nDoping1*1e6; % Doping concentration to m^-3 m2 = 0.07; % Effective electron mass coefficient m0 = 9.109382e-31; % Bare electron mass (kg) wLO=36.3; % longitudinal optic phonon energy (meV) wTO=33.3; % transverse optic phonon energy (meV) %hbar1=['\fontname{MT Extra}h\fontname{Arial}']; a0=0.529177e-10; % Bohr radius (m) einf=11.1; % high frequency dielectric constant m = m2*m0; % Effective electron mass (kg) echarge = 1.6021764e-19; % Electron charge (C) EFermi = hbar^2*(3*pi^2*nDoping)^(2/3)/(2*m); % Unit :J EFermi1 = (EFermi*1e3)/echarge; % Fermi energy in meV wLO2 = (wLO/EFermi1)^2; wTO2 = (wTO/EFermi1)^2; gamma = 0.01; %Energy broadening (GAMMA / Ef) 0.03 0.01 eye = complex(0,1); %square root of minus one kFermi = (3*pi^2*nDoping)^(1/3); %Fermi wave vector (m^-1) T = 300; % Absolute temperature 300, 10, 4.2 (K) kB = 1.3806505e-23/echarge; % Boltzmann constant (eV/K) beta = 1/(kB*T); % Inverse thermal energy (eV^-1) Energymax=300; % Electron energy maximum (meV) Energystep=Energymax/60; % Electron energy step (meV) ymax=Energymax/EFermi1; % max value of energy loss, omega/Ef 18 %max value of scattered wave vector q/kf set at ymax value xmax=(sqrt(ymax))*(1+sqrt(2));%approximately 15, 6 for n=1e17, 1e18 mu = CalculateChemicalPotential(T, nDoping, m, EFermi1/1e3); %[0.01:.005:xmax] scattered wave vector step at gamma/2 x1=[gamma/2:gamma/2:xmax]; y=[-ymax:gamma/2:ymax]+eye*gamma; % energy loss step is gamma/2 f1 = zeros(length(y),length(x1)); %************************************************************************* % Calculate the response and Im_Xe map %************************************************************************* factor = 4*pi*8.85e-12; E_q_factor=hbar^2/(2*m*echarge); Im_Xe_factor=2*echarge^2*m^2/(hbar^4*beta/echarge); Energy = y*EFermi1; %energy loss in meV for ii=1:length(x1) x = x1(ii); q = x*kFermi; %wave vector in m^-1 E_q = E_q_factor*q^2; Energy1 = (real(Energy)/1e3-E_q).^2/(4*E_q); Energy2 = (real(Energy)/1e3+E_q).^2/(4*E_q); Im_Xe = (Im_Xe_factor/(q^3))*log((1+exp(-beta*(Energy1-mu)))./(1+exp(-beta*(Energy2-mu)))); Re_Xe = -imag(hilbert(Im_Xe));%use MATLAB Hilbert transform function epsi = (y/(y-eye*gamma)*(Re_Xe+(eye*Im_Xe)))/factor+(einf*((y.^2-wLO2)./(y.^2-wTO2))); f1(:,ii)=-imag(1./epsi)./x; %loss function in the integrand end %************************************************************************* % Bose distribution %************************************************************************* Bose_mask = 1.0+(1./(exp(beta*real(y)*EFermi1/1000)-1.0)); Bose_indicator = Bose_mask'*ones(1,length(x1)); Emax_tmp = 0.1:Energystep:Energymax; spec_e = zeros(length(Emax_tmp),1); spec_eq = zeros(length(Emax_tmp),1); %************************************************************************* % Integration for electron scattering rate %************************************************************************* for ii = 1:length(Emax_tmp) Emax = Emax_tmp(ii); E_final = Emax-real(y)*EFermi1; f_fermi = 1./(exp(beta*(E_final/1000-mu))+1.0); Emax1 = Emax/EFermi1; yplot=Emax1-(x1-sqrt(Emax1)).^2; %parabola of integration yplot2=Emax1-(x1+sqrt(Emax1)).^2; %lower boundary of integration %************************************************************************* %electron-only scattering %************************************************************************* Fermi_indicator = (1.0-f_fermi')*ones(1,length(x1)); f2 = f1.*Bose_indicator.*Fermi_indicator; %************************************************************************* %quasi-particle scattering %************************************************************************* Fermi_indicatorq = (f_fermi')*ones(1,length(x1)); f2q = f1.*(Bose_indicator-Fermi_indicatorq); %************************************************************************* f2(isnan(f2))=0; %Array elements that are NaN f2q(isnan(f2))=0; %Array elements that are NaN [X1,Y1]=meshgrid(x1,y); %find values of sum indicator for each column sumindicator=Y1<=ones(length(y),1)*yplot&Y1>=ones(length(y),1)*yplot2; Integral=sum(sum(f2(sumindicator)))/(sqrt(Emax1)); Integralq=sum(sum(f2q(sumindicator)))/(sqrt(Emax1)); spec_eq(ii)=Integralq; %quasi-particle scattering rate spec_e(ii) =Integral; %electron scattering rate end spec = spec_e; %quasi-particle scattering rate energy spectrum specq = spec_eq; %electron scattering rate energy spectrum %************************************************************************* dx1 = x1(2)-x1(1); %increment in scattered wave vector X dy1 = y(2)-y(1); %increment in energy loss Y prefactor=(2.0*m2*dx1*dy1*EFermi1*0.001*echarge)/(pi*a0*kFermi*hbar); %************************************************************************* figure(1); plot(Emax_tmp,spec*prefactor/1e12,'-r'); hold on plot(Emax_tmp,specq*prefactor/1e12,'-b'); ttl3=['\rmFig7.14a, GaAs, \itn\rm_0=',num2str(nDoping1,'%4.1e'),' cm^{-3}, \itk\rm_F=',... num2str(kFermi*1e-3,'%4.1e'),' m^{-1}, \itm\rm^*_e=',num2str(m2,'%5.2f'),... '\times\itm\rm_0, \gamma=',num2str(gamma,'%5.2f'),'\times\itE\rm_F, \itT\rm=',... num2str(T,'%5.0f'),' K, \mu=',num2str(mu*1.0e3,'%5.1f'),' meV']; title(ttl3); ylabel('Inelastic electron scattering rate, 1/ \tau_{in} (ps^{-1})'); xlabel('Initial electron energy, \itE_k\rm (meV)'); axis([0,Energymax,0,30]); hold off