%Fig6_5b.m
%Essential Electron Transport for Device Physics
%Inelastic LO-phonon scattering rate at finite temperature in GaAs
clear all;
clf;

FS      = 12;	%label fontsize 18
FSN     = 12;	%number fontsize 16
LW      = 2;	%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)
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{Times}'];
a0=0.529177e-10;        % Bohr radius (m)

einf =11.1;             % high frequency dielectric constant
meff = m2*m0;           % Effective electron mass (kg)
echarge = 1.6021764e-19;            % Electron charge (C)
kscale=3.0937e+08;      %wave vector scale (m^-1)
Escale = (hbar^2)*(kscale^2)/(2*meff);% energy scale (J)
Escale1 = (Escale*1e3)/echarge;       % energy scale in meV
wLO2=(wLO/Escale1)^2;
wTO2=(wTO/Escale1)^2;
gamma=0.005;             %Energy broadening (GAMMA / Escale1) 0.01 0.005
eye=complex(0,1);        %square root of minus one

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) 100 300
Energystep=Energymax/100;           % Electron energy step (meV) /100
ymax=Energymax/Escale1;             % max value of energy loss, omega/Ef 18
%max value of scattered wave vector q/kscale set at ymax value
xmax=(sqrt(ymax))*(1+sqrt(2));%approximately 15, 6 for n=1e17, 1e18
x1=[gamma/4:gamma/4:xmax];%[0.01:.005:xmax] scattered wave vector step at gamma/4
y=[-ymax:gamma/4:ymax]+eye*gamma;% energy loss step is gamma/4
f1 = zeros(length(y),length(x1));
%*************************************************************************
% Calculate the response with phonons
%*************************************************************************
for ii=1:length(x1)
    x = x1(ii);
    epsi = einf*(1.0+((wTO2-wLO2)./(y.*(y-eye*gamma)-wTO2)));
    f1(:,ii)=-imag(1./epsi)./x; %loss function in the integrand weighted by 1/x
end
%*************************************************************************
% Bose distribution
%*************************************************************************
Bose_mask = 1.0+(1./(exp(beta*real(y)*Escale1/1000)-1.0));
Bose_indicator = Bose_mask'*ones(1,length(x1));

Emax_tmp = 0.1:Energystep:Energymax;
spec_e = zeros(length(Emax_tmp),1);
%*************************************************************************
% Integration for LO-phonon scattering rate
%*************************************************************************
%    f2 = f1.*Bose_indicator;
%*************************************************************************
f2 = f1.*Bose_indicator;
f2(isnan(f2))=0;%Array elements that are NaN
[X1,Y1]=meshgrid(x1,y);
for ii = 1:length(Emax_tmp)
    Emax1 = Emax_tmp(ii)/Escale1;

    yplot =Emax1-(x1-sqrt(Emax1)).^2;  %parabola of integration
    yplot2=Emax1-(x1+sqrt(Emax1)).^2; %lower boundary of integration

%quasi-particle scattering    
    sumindicator=Y1<=ones(length(y),1)*yplot&Y1>=ones(length(y),1)*yplot2; %find values of sum indicator for each column
    Integral=sum(sum(f2(sumindicator)))/(sqrt(Emax1));  
    spec_e(ii)=Integral;
end
spec = spec_e;

%*************************************************************************
figure(1);
surf(x1*kscale,real(y)*Escale1,log10(abs(f1)));
                az = 0;           %set up the angle of view
                el = 90;          %looking from directly above
                view(az, el);     %set the view
                shading interp;   %do color interpolation
       colorbar;
       axis([0 xmax*kscale -Energymax Energymax]);
       hold on;
       plot(x1*kscale, yplot*Escale1,'r','LineWidth',2);
       plot(x1*kscale, yplot2*Escale1,'r','LineWidth',2);
       hold off;
       xlabel('Wave vector, \itq\rm (m^{-1})');
       ylabel(['Energy transfer, ',hbar1,'\omega (meV)']);
ttl1=['\rmGaAs log_{10}(\itS\rm(\itq\rm, \omega )/\itq\rm), \itk\rm_{scale}=',num2str(kscale,...
    '%4.1e'),' m^{-1}, \itm\rm^*_e=',num2str(m2,'%5.2f'),...
    '\times\itm\rm_0, \gamma=',num2str(gamma*Escale1,'%5.2f'),...
    ' meV, \itT\rm=',num2str(T,'%5.1f'),' K'];
title(ttl1);
%*************************************************************************
figure(2);
surf(x1*kscale,real(y)*Escale1,log10(abs(f2)));
                az = 0;           %set up the angle of view
                el = 90;          %looking from directly above
                view(az, el);     %set the view
                shading interp;   %do color interpolation
       colorbar;
       axis([0 xmax*kscale -Energymax Energymax]);
       hold on;
       plot(x1*kscale, yplot*Escale1,'w','LineWidth',LW);
       plot(x1*kscale, yplot2*Escale1,'w','LineWidth',LW);
       hold off;
       xlabel('Wave vector, \itq\rm (m^{-1})');
       ylabel(['Energy transfer, ',hbar1,'\omega (meV)']);
ttl2=['\rmGaAs log_{10}(\itS\rm(\itq\rm, \omega )\itg\rm/\itq\rm), \itk\rm_{scale}=',num2str(kscale,...
    '%4.1e'),' m^{-1}, \itm\rm^*_e=',num2str(m2,'%5.2f'),...
    '\times\itm\rm_0, \gamma=',num2str(gamma*Escale1,'%5.2f'),...
    ' meV, \itT\rm=',num2str(T,'%5.1f'),' K'];
title(ttl2);
%*************************************************************************
dx1 = x1(2)-x1(1);
dy1 = y(2)-y(1);
prefactor=(2.0*m2*dx1*dy1*Escale1*0.001*echarge)/(pi*a0*kscale*hbar);
%*************************************************************************
figure(3);
plot(Emax_tmp,spec*prefactor/1e12);
ttl3=['\rmFig6.5b, GaAs, \itk\rm_{scale}=',num2str(kscale,...
    '%4.1e'),' m^{-1}, \itm\rm^*_e=',num2str(m2,'%5.2f'),...
    '\times\itm\rm_0, \gamma=',num2str(gamma*Escale1,'%5.2f'),...
    ' meV, \itT\rm=',num2str(T,'%5.1f'),' K'];
title(ttl3);
ylabel('Inelastic electron scattering rate, 1/\tau_{in} (ps^{-1})');
xlabel('Initial electron energy, \itE_k\rm (meV)');
grid on;
axis([0,Energymax,0,10]);
max(spec*prefactor/1e12)
%*************************************************************************

