%Fig6_3a.m
%Essential Electron Transport for Device Physics
%calculates optical phonon scattering rates
%eg GaAs meff=0.067 wLO=36.3, wTO=33.3, epsi_inf=11.1
%InAs meff=0.021 wLO=30.2, wTO=27.1, epsi_inf=11.8
%InP  meff=0.077 wLO=42.7, wTO=38.2, epsi_inf=9.52 
%
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{Times}'];
mfpl_txt=['\fontname{MT Extra}l\fontname{Times}'];
%**************************************************************************
m0=9.1095e-31;                  %bare electron mass [kg]
meff=0.07;                      %effective electron mass in conduction band
wLO=36.3;                       %longitudinal optic phonon energy [meV]
wTO=33.3;                       %transverse optic phonon energy [meV]
epsi_inf=11.1;                  %high frequency dielectric constant
epsi_0=epsi_inf*((wLO/wTO)^2);  %low frequency dielectric constant
hbar=1.05459e-34;               %Plank constant [J s]
e0=8.8541878e-12;               %Permittivity of free space [F m^-1]
echarge=1.60219e-19;            %electron charge [C]
kB=1.3806503e-23;               %Boltzmann constant [J K^-1]
%**************************************************************************
%	epsilon zero, epsilon inf., electron mass, phonon frequency
%**************************************************************************

    emax=300;       %maximum electron energy [meV]
    npoints=300;    %number of points in plot
    
	tmin =0.0;      %minimum temperature [K]
    tstep=100.0;    %temperature step [K]
	tmax =500;      %max temperature [K]
%
%   Set up plot coordinates
%
	xmin=0.0;       %minimum energy in plot
	xmax=emax;      %maximum energy in plot [meV]
	ymin=0.0;       %minimum scattering rate in plot
	ymax=1.5e13;    %maximum scattering rate in plot
%	
% setup prefactor in scattering rate calaculation
%
    g0=m0*meff*(echarge^2)*(echarge*wLO/1000)/(4*pi*e0*(hbar^2));
    g0=g0/(sqrt(2*m0*meff*echarge));
    g0=g0*((1.0/epsi_inf)-(1.0/epsi_0));    
%******************** loop *******************   
eEmeV=linspace(0.0,emax,npoints);           %electron energy [meV]
eE=eEmeV/1000;                              %electron energy [eV]

for Ttemp=tmin:tstep:tmax                   %main loop in temperature
beta=1/(kB*Ttemp);                          %inverse thermal energy in J
nq=1.0/(exp(beta*wLO*echarge/1000)-1.0);    %Bose thermal occupation factor
%  tmp1=abs((sqrt(eEmeV+wLO)+sqrt(eEmeV))./(sqrt(eEmeV+wLO)-sqrt(eEmeV)));
%  tmp2=abs((sqrt(eEmeV)+sqrt(eEmeV-wLO))./(sqrt(eEmeV)-sqrt(eEmeV-wLO))); 
%  g1=nq*log(tmp1);%phonon absorption
%  g2=(nq+1.0)*log(tmp2);%phonon emission
tmp1=2*asinh(sqrt(eEmeV/wLO));
tmp2=2*asinh(sqrt((eEmeV-wLO)/wLO));
tmp2=real(tmp2);%remove terms when eEmeV<wLO
g1=nq*tmp1;                 %phonon absorption
g2=(nq+1.0)*tmp2;           %phonon emission    
%     
Gamma =g0*(g1+g2)./sqrt(eE);%total scattering rate
Gamma1=g0*g1./sqrt(eE);     %scattering rate for absorption of phonons only     
%********** plot scattering rate figure ***********************************    
figure(1);   
plot(eEmeV,Gamma1*1e-12,':r','LineWidth',LW); %absorption scattering rate
grid on
hold on
plot(eEmeV,Gamma*1e-12,'-b','LineWidth',LW);  %total scattering rate
end
max(Gamma*1e-12)
%**************************************************************************    
ttl1=['\rmFig6.3a, GaAs LO-phonon, \itm\rm^*_e=',num2str(meff),...
'\times\itm\rm_0, ',hbar_txt,'\omega_{LO}=',num2str(wLO),...
' meV, T_{min}=',num2str(tmin),' K, \itT\rm_{step}=',num2str(tstep),' K'];    
xlabel('Energy, \itE_k\rm (meV)');
ylabel('Scattering rate, 1/ \tau_{LO} (ps^{-1})');
axis([xmin xmax ymin ymax*1e-12]);
title(ttl1);
yscale=ylim;
y1=yscale(2)/1.1;
y2=y1*0.1;
text(xmax*0.72,y1,['T_{max} = ',num2str(tmax),' K']);
text(xmax*0.72,y2,['T_{min} = ',num2str(tmin),' K']);
   
%velocity [m s^-1] of electron energy eE [eV]
vel=sqrt(2*eE*echarge/(m0*meff));
for Ttemp=tmin:tstep:tmax                   %main loop in temperature
beta=1/(kB*Ttemp);                          %inverse thermal energy in J
nq=1.0/(exp(beta*wLO*echarge/1000)-1.0);    %Bose thermal occupation factor
    
%  tmp1=abs((sqrt(eEmeV+wLO)+sqrt(eEmeV))./(sqrt(eEmeV+wLO)-sqrt(eEmeV)));
%  tmp2=abs((sqrt(eEmeV)+sqrt(eEmeV-wLO))./(sqrt(eEmeV)-sqrt(eEmeV-wLO))); 
%  g1=nq*log(tmp1);%phonon absorption
%  g2=(nq+1.0)*log(tmp2);%phonon emission
tmp1=2*asinh(sqrt(eEmeV/wLO));
tmp2=2*asinh(sqrt((eEmeV-wLO)/wLO));    
tmp2=real(tmp2);%remove terms when eEmeV<wLO
g1=nq*tmp1;                 %phonon absorption
g2=(nq+1.0)*tmp2;           %phonon emission    
%     
Gamma =g0*(g1+g2)./sqrt(eE);%total scattering rate
mfp = vel./Gamma; %calculate mean-free-path [m]

%********** plot mean-free-path figure ************************************    
figure(2);   
grid on
hold on
plot(eEmeV,mfp*1e9,'-b','LineWidth',LW);  %mean-free-path in nm
end
%**************************************************************************    
hold off
ttl1=['\rmGaAs LO-phonon mean-free-path, \itm\rm^*_e=',num2str(meff),...
'\times\itm\rm_0, ',hbar_txt,'\omega_{LO}=',num2str(wLO),...
' meV, \itT\rm_{min}=',num2str(tmin),' K, \itT\rm_{step}=',num2str(tstep),' K'];
ttl2=['\rmMean-free-path, ',mfpl_txt,' (nm)'];
xlabel('Energy, \itE_k\rm (meV)');
ylabel(ttl2);
axis([xmin xmax ymin 500]);
title(ttl1);