% Calculate electron distribution of GaSb with nonparabolic bands 
% in different conduction band minima 
% as a function of electron concentration at different temperatures.
% Prepared by Stefan Zollner, New Mexico State University, Las Cruces, NM.
% Date: 14 June 2024
% Required package: Symbolic Math Toolbox Add-On, for polylog() function.
% This script takes about ten minutes to run on a Windows 11 laptop.
%
% References:
% S. Zollner, S. Gopalan, and M. Cardona, 
% The temperature dependence of the band gaps in InP, InAs, InSb, and GaSb, 
% Solid-State Commun. 77, 485-488 (1991). https://doi.org/10.1016/0038-1098(91)90725-B
%
% J. Camassel and D. Auvergne, 
% Temperature dependence of the fundamental edge of germanium and 
% zinc-blende-type semiconductors,
% Phys. Rev. B 12, 3258 (1975). https://doi.org/10.1103/PhysRevB.12.3258
%
% S. C. Chen and Y. K. Su,
% Photoluminescence study of gallium antimonide grown by liquid‐phase epitaxy,
% J. Appl. Phys. 66, 350–353 (1989). https://doi.org/10.1063/1.343880
%
% Landolt-Boernstein, Numerical Data and Functional Relationships in
% Science and Technology, Volume 22a and 17a. 
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Physical constants (energy units in eV unless noted) from CODATA at NIST
hbar=6.582119569E-16; % eV*s
e=1.602176634E-19; % A*s
m0=9.1093837015E-31; % kg
% eps0=8.8541878128E-12; % F/m
kB=0.08617333262E-3; % eV/K Boltzmann constant
% kBJoule=kB*e;
% c=299792458; % m/s speed of light in vacuum
%
% h=hbar*2*pi; % eV*s
% hJoule=h*e; % J*s
hbarJoule=hbar*e;
% hbarc=hbar*c*100; % eV*cm
% hbar2m0=hbar*hbarJoule/m0*1E20; % hbar^2/m0 in Angstrom^2
% Rydberg=e^4*m0/(8*eps0^2*hJoule^2)/e; % eV, hydrogen ground state energy
% Hartree=2*Rydberg; % eV, atomic unit for energy 
% aBohr=(hJoule^2*eps0)/(pi*e^2*m0)*1E10; % Angstrom
%Prefactor e^2/(eps0*m0) for plasma frequency
% PlasmaFactor=e^2/(eps0*m0)*1E6; % cm-3/s^2
% DebyeFactor=100*e^2/(eps0*kBJoule); % Kcm
% amplitude prefactor for optical interband transitions
% Afactor=4*sqrt(Rydberg)/3;
% prefactor for the electron and hole concentrations in cm^{-3}
Nfactor=((2*m0*kB)/(pi*hbar*hbarJoule))^1.5*1E-6/4;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Material parameters for GaSb (compare Ioffe Institute web site)
a=6.09593; % lattice constant of GaSb in Angstrom (Landolt-Boernstein vol. 22a)
%
% ThetaD=450; % Parameter for Debye temperature (266 K in GaSb, compare LaBo vol. 22a)
ThetaTA=0.00703/kB; % TA phonon energy at the X point (in eV)
cLA=3.97E3; % longitudinal speed of sound in m/s
% calculated as sqrt(c11/rho), c11=88.34, rho=5.6137, see LaBo vol. 22a.
ThetaLA=hbar*cLA/kB*2*pi/a/1E-10*nthroot(3/pi,3);
omegaR=0.0284;         % first order Raman energy in eV for InSb (Ioffe)
ThetaO=0.9*omegaR/kB; % See Eq. (10) in Roucka
ATA=-2.53;
ALA=5.87;
AO=ALA;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
me=0.041; % effective mass at Gamma at 0 K (LaBo 22a)
mLt=0.11; % transverse mass at L (LaBo 22a)
mLl=0.95; % longitudinal mass at L (LaBo 22a)
mL=nthroot(mLt*mLt*mLl,3); % DOS mass at L 
mXt=0.22; % transverse mass at X (LaBo 22a)
mXl=1.51; % longitudinal mass at X (LaBo 22a
mX=nthroot(mXt*mXt*mXl,3); % DOS mass at X
E0=0.811; % direct gap of GaSb at low T in eV (Landolt-Boernstein vol. 22a)
DeltaGL=0.089; % difference between Gamma and L valley energies in eV (LaBo 22a)
DeltaGX=0.33 ; % difference between Gamma and X valley energies in eV (LaBo 22a)
E1=2.184; % E1 energy in eV, from Zollner Ph.D. thesis, Table 9.1
Delta1=0.434; % spin-orbit splitting at L in eV, from Zollner Ph.D. thesis, Table 9.1
alphaep=(1-me)^2; % nonparabolicity parameter
alphae=alphaep/E0;
DeltaL=1/(1/E1+1/(E1+Delta1)); % nonparabolicity parameter at L from Menendez PL 2020
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% E0T=@(T)E0;     % keep E0 constant at value given above
E0T1=@(T)0.810-3.78E-4*T.^2./(T+94); 
% from Camassel and Auvergne, PRB 12, 3258 (1975).
% as digitized in Casey and Panish, Heterostructure Lasers, Part B, page 9.
% also cited by Chen and Su, JAP 66, 350 (1989).
E0T2=@(T)0.812-4.2E-4*T.^2./(T+140);
% from Chen and Su, JAP 66, 350 (1989).
% Bose-Einstein parameters for temperature dependence of experimental band gap
EBE0=0.831; % eV
aBE0=0.02; % eV
OmegaE0=0.01; % eV
% Full temperature dependence of the experimental direct band gap using Bose-Einstein expression.
% This gap E0BE includes thermal expansion and electron-phonon contributions. 
E0T=@(T)(EBE0 - aBE0* (1 + 2./(exp(OmegaE0./(kB*T))-1)) ); 
% The parameters were fitted to Camassel and Auvergne
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Fermi=@(m,x)-polylog(m+1,-exp(x)); % Fermi-Dirac integrals expressed as polylogs.
% Fermi=@(m,x)exp(x); classical Maxwell-Boltzmann approximation
options=optimset('TolX',1E-4);
%
Tstart=10;
Tstep=10;
Nstep=30;
for i=2:Nstep+1
    Ti(i)=Tstart+Tstep*(i-2);
end 
%
GaSbTemp=figure('Name','Temperature dependence of the band gap of GaSb', ...
             'Color','white');
pl=plot(Ti,E0T1(Ti),'r--',Ti,E0T2(Ti),'b:',Ti,E0T(Ti),'k-');
title('Temperature dependence of the band gap of GaSb');
xlabel('Temperature T (K)'); ylabel('Band Gap E_0 (eV)');
xlim([0,300]); % ylim([-1,1]);
pl(1).LineWidth=2; pl(2).LineWidth=2; pl(3).LineWidth=2; 
legend('Camassel & Auvergne 1975','Chen & Su 1989','Bose-Einstein fit', ...
      'FontSize',14, 'Location','southwest'); 
legend('boxoff');
% close(GaSbTemp); 
% return;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
muvsT=figure('Name','Chemical potential versus temperature', ...
             'Color','white');
%
jstart=10; jend=19; jstep=1;
nGfrac=zeros(Nstep+1,jend); 
nLfrac=zeros(Nstep+1,jend); 
nXfrac=zeros(Nstep+1,jend); 
%
for j=jstart:jstep:jend
    nlaser=10^j;
    fprintf('n=%5.0E cm^{-3}\n',nlaser);
    Ti=zeros(1,Nstep+1); 
%    muroot=zeros(1,Nstep+1); murootnp=zeros(1,Nstep+1); 
    murootne=zeros(1,Nstep+1); murootnenp=zeros(1,Nstep+1);
%    
    for i=2:Nstep+1
        Ti(i)=Tstart+Tstep*(i-2);
        T=Ti(i);   
        if i==2 
            muguess=E0T(T); muguessnp=muguess;
        else 
%             muguess=muroot(i-1); muguessnp=murootnp(i-1);
            muguess=murootne(i-1); muguessnp=murootnenp(i-1);
        end
        neG=@(x)Nfactor*(me*T)^1.5*Fermi(0.5,(x-E0T(T))/(kB*T));
        neGnp=@(x)neG(x)+Nfactor*(me*T)^1.5*15/4*alphae*kB*T*Fermi(1.5,(x-E0T(T))/(kB*T));
        nL=@(x)4*Nfactor*(mL*T)^1.5*Fermi(0.5,(x-E0T(T)-DeltaGL)/(kB*T));
        nLnp=@(x)nL(x)+4*Nfactor*(mL*T)^1.5*kB*T/DeltaL*Fermi(1.5,(x-E0T(T)-DeltaGL)/(kB*T));
% Multiply by four (number of L-valleys)
        nX=@(x)3*Nfactor*(mX*T)^1.5*Fermi(0.5,(x-E0T(T)-DeltaGX)/(kB*T));
% Multiply by three (number of X-valleys). This would six for diamond-type materials.
        ne=@(x)neG(x)+nL(x)+nX(x);
        nenp=@(x)neGnp(x)+nLnp(x)+nX(x);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%        diff=@(x)neG(x)-nlaser;
%        diffnp=@(x)neGnp(x)-nlaser;
        diffne=@(x)ne(x)-nlaser;
        diffnenp=@(x)nenp(x)-nlaser;
%        muroot(i)=fzero(diff,muguess,options);
%        murootnp(i)=fzero(diffnp,muguessnp,options);
%        fprintf('T=%3.0f K, muG=%5.3f eV, muGnp=%5.3f eV\n', ...
%                T,muroot(i),murootnp(i));
        murootne(i)=fzero(diffne,muguess,options);
        murootnenp(i)=fzero(diffnenp,muguessnp,options);
        fprintf('T=%3.0f K, mu=%5.3f eV, munp=%5.3f eV\n', ...
                T,murootne(i),murootnenp(i));
        nGfrac(i,j)=neGnp(murootnenp(i))/nlaser;
        nLfrac(i,j)=nLnp(murootnenp(i))/nlaser;
        nXfrac(i,j)=nX(murootnenp(i))/nlaser;
        fprintf('T=%3.0f K, nGfrac=%5.3f, nLfrac=%5.3f, nXfrac=%5.3f\n', ...
                T,nGfrac(i,j),nLfrac(i,j),nXfrac(i,j));
%
    end
%
    Ti(1)=1; % initial temperature in Kelvin with special processing
    T=1;
% This section is for Gamma-valley only (ignore carriers near L and X)
%    if j<16 % use classical Maxwell-Boltzmann statistics
%        neG=@(x)Nfactor*(me*T)^1.5*exp((x-E0T(T))/(kB*T));
%        neGnp=@(x)neG(x)+Nfactor*(me*T)^1.5*15/4*alphae*kB*T*exp((x-E0T(T))/(kB*T));
%        diff=@(x)neG(x)-nlaser;
%        diffnp=@(x)neGnp(x)-nlaser;
%        muroot(1)=fzero(diff,E0T(T)+0.01,options);
%        murootnp(1)=fzero(diffnp,E0T(T)+0.01,options);
%    else
%        muroot(1)=(nlaser/Nfactor/(me*T)^1.5*3*sqrt(pi)/4)^(2/3)*kB*T+E0T(T);
%        murootnp(1)=murootnp(2); % here we cheat
%    end
%    fprintf('T=%3.0f K, muG=%5.3f eV, muGnp=%5.3f eV\n', ...
%            T,muroot(1),murootnp(1));
%    murootne(1)=murootne(2);        % here we cheat and take the value
%    murootnenp(1)=murootnenp(2);    % at the next highest temperature
% This section includes the Gamma, L, and X-valleys
    neG=@(x)Nfactor*(me*T)^1.5*Fermi(0.5,(x-E0T(T))/(kB*T));
    neGnp=@(x)neG(x)+Nfactor*(me*T)^1.5*15/4*alphae*kB*T*Fermi(1.5,(x-E0T(T))/(kB*T));
    nL=@(x)4*Nfactor*(mL*T)^1.5*Fermi(0.5,(x-E0T(T)-DeltaGL)/(kB*T));
    nLnp=@(x)nL(x)+4*Nfactor*(mL*T)^1.5*kB*T/DeltaL*Fermi(1.5,(x-E0T(T)-DeltaGL)/(kB*T));
% Multiply by four (number of L-valleys)
    nX=@(x)3*Nfactor*(mX*T)^1.5*Fermi(0.5,(x-E0T(T)-DeltaGX)/(kB*T));
% Multiply by three (number of X-valleys). This would six for diamond-type materials.
    ne=@(x)neG(x)+nL(x)+nX(x);
    nenp=@(x)neGnp(x)+nL(x)+nX(x);
    diffne=@(x)ne(x)-nlaser;
    diffnenp=@(x)nenp(x)-nlaser;
%
    if j<=17
       murootne(1)=fzero(diffne,E0T(T),options);
       murootnenp(1)=fzero(diffnenp,E0T(T),options);
       nGfrac(1,j)=neGnp(murootnenp(1))/nlaser;
       nLfrac(1,j)=nLnp(murootnenp(1))/nlaser;
       nXfrac(1,j)=nX(murootnenp(1))/nlaser;
    else
       murootne(1)=murootne(2);         % here we cheat
       murootnenp(1)=murootnenp(2);     % take the results from the next temperature
       nGfrac(1,j)=nGfrac(2,j);
       nLfrac(1,j)=nLfrac(2,j);
       nXfrac(1,j)=nXfrac(2,j);
    end
    fprintf('T=%3.0f K, muG=%5.3f eV, muGnp=%5.3f eV\n', ...
            T,murootne(1),murootnenp(1));

    fprintf('T=%3.0f K, nGfrac=%5.3f, nLfrac=%5.3f, nXfrac=%5.3f\n', ...
            T,nGfrac(1,j),nLfrac(1,j),nXfrac(1,j));
%
    pl=plot(Ti,murootne,'r-',Ti,murootnenp,'g-');
    pl(1).LineWidth=2; pl(2).LineWidth=2; 
    hold on;
end
%
hold off;
% title('GaSb: Fixed band gap E_0='+string(sprintf('%5.3f',E0))+' eV, constant masses');
title('GaSb: Bose-Einstein band gap, constant masses');
xlabel('Temperature T (K)'); ylabel('chemical potential (eV)');
xlim([0,300]); % ylim([-1,1]);
legend('parabolic bands (\Gamma+L+X)','nonparabolic bands (\Gamma+L+X)', ...
      'FontSize',14, 'Location','southwest'); 
% legend('parabolic bands (\Gamma only)','nonparabolic bands (\Gamma only)', ...
%       'Location','southwest'); 
% This text is for the figure with Gamma only (ignore L- and X-valleys)
% text(0.8,0.05,string(sprintf('%5.0E',10^jstart)), ...
%      'FontSize',14,'Units','normalized','Color','green');
% text(0.8,0.75,string(sprintf('%5.0E',10^jend)), ...
%      'FontSize',14,'Units','normalized','Color','green');
% text(0.8,0.88,string(sprintf('%5.0E',10^jend)), ...
%      'FontSize',14,'Units','normalized','Color','red');
% text(0.05,0.64,string(sprintf('%5.0E',10^18)), ...
%      'FontSize',14,'Units','normalized','Color','red');
text(0.85,0.05,string(sprintf('%5.0E',10^jstart)), ...
     'FontSize',14,'Units','normalized','Color','green');
text(0.85,0.93,string(sprintf('%5.0E',10^jend)), ...
     'FontSize',14,'Units','normalized','Color','green');
legend('boxoff');
% 
nvsT=figure('Name','Electron populations versus temperature', ...
             'Color','white');
%
for j=16:jstep:jend
    pl=plot(Ti,nGfrac(:,j),'r-',Ti,nLfrac(:,j),'g-',Ti,nXfrac(:,j),'b-');
    pl(1).LineWidth=2; pl(2).LineWidth=2; pl(3).LineWidth=2;
    hold on;
end
hold off;
% title('GaSb: Fixed band gap E_0='+string(sprintf('%5.3f',E0))+' eV, constant masses');
title('GaSb: Bose-Einstein band gap, constant masses');
xlabel('Temperature T (K)'); ylabel('relative electron population');
xlim([0,300]); ylim([0,1]);
% legend('\Gamma','L', 'X','FontSize',14,'Position',[0.2 0.5 0.2 0.2]); 
legend('\Gamma','L', 'X','FontSize',14,'Location','west'); 
legend('boxoff');
text(250,0.9,'10^{19}','FontSize',14,'Color','green');
text(250,0.7,'10^{18}','FontSize',14,'Color','green');
text(180,0.3,'10^{17}','FontSize',14,'Color','green');
text(220,0.28,'10^{16}','FontSize',14,'Color','green');
text(30,0.22,'10^{19}','FontSize',14,'Color','red');
text(120,0.6,'10^{18}','FontSize',14,'Color','red');
text(170,0.7,'10^{17}','FontSize',14,'Color','red');
text(215,0.75,'10^{16}','FontSize',14,'Color','red');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
return;
