%Fig6_3b.m
%Essential Electron Transport for Device Physics
%calculates optical phonon scattering rates for fixed electron energy
%as function of temperature, T
%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}'];
%**************************************************************************
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
%**************************************************************************

	tmin=0.0;       %minimum temperature in K
    tstep=1.0;      %temperature step in K
	tmax=500;       %max temperature in K
%
%   Set up plot coordinates
%
	xmin=0.0;       %minimum energy in plot
	xmax=tmax;      %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=[wLO/2,100,200,300];                  %electron energy in meV
eE=eEmeV/1000;                              %electron energy in eV
Ttemp=linspace(tmin,tmax,500);
%
for n=1:length(eEmeV);                   %main loop in energy
beta=1./(kB*Ttemp);                      %inverse thermal energy in J
nq=1.0./(exp(beta*wLO*echarge/1000)-1.0);%Bose thermal occupation factor
    
%keE=sqrt(2*m0*ms*eE)/hb;%electron wave vector in m^-1
    
%  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(n)/wLO));
tmp2=2*asinh(sqrt((eEmeV(n)-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(n));%total scattering rate
Gamma1=g0*g1./sqrt(eE(n));     %scattering rate for absorption of phonons only     
%********** plot scattering rate figure ***********************************    
figure(1);   
plot(Ttemp,Gamma1*1e-12,':r','LineWidth',LW); %absorption scattering rate
grid on
hold on
plot(Ttemp,Gamma*1e-12,'b','LineWidth',LW);  %total scattering rate
end
%**************************************************************************    
ttl1=['\rmFig6.3b, GaAs LO-phonon, \itm\rm^*_e = ',num2str(meff),...
'\times\itm\rm_0, ',hbar_txt,'\omega_{LO} = ',num2str(wLO),' meV'];    
xlabel('Temperature, \itT\rm (K)');
ylabel('Scattering rate, 1/ \tau_{LO} (ps^{-1})');
axis([xmin xmax ymin ymax*1e-12]);
title(ttl1);
yscale=ylim;
y4=yscale(2)/1.1;
y3=y4*0.9;
y2=y4*0.8;
y1=y4*0.7;
text(xmax*0.05,y4,['E_{4} = ',num2str(eEmeV(4)),' meV']);
text(xmax*0.05,y3,['E_{3} = ',num2str(eEmeV(3)),' meV']);
text(xmax*0.05,y2,['E_{2} = ',num2str(eEmeV(2)),' meV']);
text(xmax*0.05,y1,['E_{1} = ',num2str(eEmeV(1)),' meV']);