%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
