% Calculate thermal expansion coefficient of InSb.
% Also calculate redshift of the band gap due to thermal expansion.
% Prepared by Stefan Zollner, New Mexico State University, Las Cruces, NM.
% Date: 29 December 2023
% Required package: None.
%
% References:
% Stefan Zollner, Carlos A. Armenta, Sonam Yadav, and Jose Menendez, 
% Conduction band nonparabolicity, chemical potential, abd carrier
% concentration of intrinsic InSb as a function of temperature,
% J. Vac. Soc. Technol. MIOMD-XI special collection (submitted 2024). 
%
% R. Roucka, Y.-Y. Fang, J. Kouvetakis, A. V. G. Chizmeshya, and J. Menendez,
% Thermal expansivity of Ge1-ySny alloys, Phys. Rev. B 81, 245214 (2010).
% S. I. Novikova, Thermal expansion, in Semiconductors and Semimetals, vol. 2, 
% edited by R. K. Willardson and A. C. Beer, pp. 33-38 (1966).
% D. F. Gibbons, Thermal expansion of some crystals with the diamond
% structure, Phys. Rev. 112, 136 (1958). 
% P. W. Sparks and C. A. Swenson, Thermal expansion from 2 to 40K of Ge,
% Si, and four III/V compounds, Phys. Rev. 163, 779 (1967).
% X. Cai and J. Wei, Temperature dependence of the thermal properties of
% InSb materials used in data storage, J. Appl. Phys. 114, 083507 (2013). 
% Stefan Zollner, Sudha Gopalan, and Manuel Cardona, The temperature dependence of the 
% band gaps in InP, InAs, InSb, and GaSb, Solid State Commun. 77, 485-488 (1991). 
% Stefan Zollner, Miquel Garriga, Josef Humlicek, Sudha Gopalan, and Manuel Cardona,
% Temperature dependence of the dielectric function and the interband critical-point
% parameters of GaSb, Phys. Rev. B 43, 4349-4360 (1991). 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% experimental linear expansion coefficients for InSb 
% as digitized in Stuttgart in 1991, taken from S.I. Novikova, in 
% Semiconductors and Semimetals, edited by R.K. Willardson and A.C. Beer,
% vol. 2, Academic Press, New York, 1966.
% in units of 10^{-6} K^{-1}
alphaexp=[ ...
  20.000     -1.6; ...
  30.000     -1.5; ...
  40.000     -1.2; ...
  50.000     -0.6; ...
  60.000     0.2; ...
  70.000      1.0000; ...
  80.000      1.6000; ...
  90.000      2.3000; ...
  100.00      2.8000; ...
  110.00      3.000; ...
  120.00      3.2000; ...
  130.00      3.3000; ...
  140.00      3.8000; ...
  150.00      4.0000; ...
  160.00      4.2000; ...
  170.00      4.3000; ...
  180.00      4.4000; ...
  190.00      4.4500; ...
  200.00      4.5000; ...
  210.00      4.5500; ...
  220.00      4.6000; ...
  230.00      4.6500; ...
  240.00      4.7500; ...
  250.00      4.8500; ...
  260.00      5.0000; ...
  270.00      5.0500; ...
  280.00      5.1000; ...
  290.00      5.1500; ...
  300.00      5.2000];
% Data from Xiaolin Cai and Jingsong Wei, J. Appl. Phys. 114, 083507 (2013)
alphaCai=[ ...
    300 4.8; ...
    370 5.0; ...
    460 5.5; ...
    570 5.8; ...
    660 5.6; ...
    770 5.5];
%
alphaGibbons=[ ...
20	-0.1;   ...
30	-1.72;  ...
40	-0.82;  ...
50	-0.33;  ...
60  0.28;   ...
70	0.89;   ...
80	1.5;    ...
90	2.18;   ...
100	2.76;   ...
110	2.22;   ...
120	2.48;   ...
130	3.67;   ...
140	3.83;   ...
150	3.96;   ...
160	4.08;   ...
170	4.17;   ...
180	4.27;   ...
190	4.35;   ...
200	4.43;   ...
210	4.51;   ...
220	4.58;   ...
230	4.64;   ...
240	4.71;   ...
250	4.78;   ...
260	4.83;   ...
270	4.89;   ...
280	4.95;   ...
290	5.00;   ...
300	5.04    ];
% 
alphaSparks=[ ...
6   0;...
10  -0.25;...
15  -0.8;...
20  -1.2;...
25  -1.4;...
30  -1.38;...
35  -1.2;];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a=6.479; % lattice constant of InSb in Angstrom (Ioffe Institute)
ThetaD=450; % Parameter for Debye temperature (205 K in InSb, compare LaBo vol. 22a)
ThetaTA=0.00463/kB;
cLA=3.4E3; % longitudinal speed of sound in m/s
ThetaLA=hbar*cLA/kB*2*pi/a/1E-10*nthroot(3/pi,3);
omegaR=0.023;         % first order Raman energy in eV for InSb (Ioffe)
ThetaO=0.9*omegaR/kB; % See Eq. (10) in Roucka
ALA=0.055;
ATA=-3;
AO=7;
amp=17.5; % adjustable parameter to fit the magnitude of the thermal expansion coeff. 
fun=@(x)x.^4.*exp(x)./(exp(x)-1).^2; % function under the Debye integral
%
Tvar=zeros(1,80); alphavar=zeros(1,80); alphavar8=zeros(1,80); 
alphaint=zeros(1,80); alphaint8=zeros(1,80);  
Tstart=10;
Tstep=10;
alphaint(1)=0; 
for i=1:1:80
    Ti=Tstart+(i-1)*Tstep;
    Tvar(i)=Ti; 
% Debye integral solved numerically
    alphavar(i)=amp*(Ti/ThetaD)^3*integral(fun,0,ThetaD/Ti); 
    alphavar8(i)=ATA*(ThetaTA/Ti)^2*exp(ThetaTA/Ti)/(exp(ThetaTA/Ti)-1)^2+ ...
                 ALA*(Ti/ThetaTA)^3*integral(fun,0,ThetaLA/Ti)+ ...
                 AO*(ThetaO/Ti)^2*exp(ThetaO/Ti)/(exp(ThetaO/Ti)-1)^2;
% Equation (6) from Roucka, PRB 81, 245214 (2000).
% calculate integral over alpha(t)dT from temperature 0 to T.    
    if i~=1 
        alphaint(i)=alphaint(i-1)+alphavar(i)*Tstep; 
        alphaint8(i)=alphaint8(i-1)+alphavar8(i)*Tstep;
    end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f=figure('Name','Linear thermal expansion coefficient of InSb', ...
    'Color','white');
pl=plot(Tvar,alphavar,':',Tvar,alphavar8,'-');
pl(1).LineWidth=2; pl(2).LineWidth=2; 
hold on;
pl=plot(alphaexp(1:29,1),alphaexp(1:29,2),'+r', ...
        alphaCai(1:6,1),alphaCai(1:6,2),'*k', ...
        alphaGibbons(1:28,1),alphaGibbons(1:28,2),'og', ...
        alphaSparks(1:7,1),alphaSparks(1:7,2),'xc');
hold off;
ax=gca;
ax.XAxis.FontSize=12;
ax.YAxis.FontSize=12;
xlabel('Temperature T (K)',FontSize=12);
% ylabel('linear thermal expansion coefficient \alpha (10^{-6} K^{-1})'); 
ylabel('\alpha (10^{-6} K^{-1})'); 
pl(1).LineWidth=2; 
legend('Roucka Eq. (6)','Roucka Eq. (8)', ...
       'Novikova','Cai','Gibbons','Sparks', ...
       'Location','southeast','FontSize',14); 
legend('boxoff');
text(0.2,0.3,'A='+string(amp)+'\times{}10^{-6} K^{-1}', ...
    'FontSize',14,'Units','normalized');
text(0.2,0.2,'\Theta_D='+string(ThetaD)+' K', ...
    'FontSize',14,'Units','normalized');
text(0.5,0.8,'A_{\rm TA}='+string(ATA)+'\times{}10^{-6} K^{-1}', ...
    'FontSize',14,'Units','normalized');
text(0.5,0.7,'A_{\rm LA}='+string(ALA)+'\times{}10^{-6} K^{-1}', ...
    'FontSize',14,'Units','normalized');
text(0.5,0.6,'A_O='+string(AO)+'\times{}10^{-6} K^{-1}', ...
    'FontSize',14,'Units','normalized');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate redshift of the direct band gap of InSb due to thermal expansion. Compare: 
% Stefan Zollner, Sudha Gopalan, and Manuel Cardona, The temperature dependence of the 
% band gaps in InP, InAs, InSb, and GaSb, Solid State Commun. 77, 485-488 (1991). 
% E0=0.237; % eV, from Lawaetz 1971, at 0 K
B=46; % bulk modulus in GPa, LaBoe Vol. 17a, page 317
dEdp=0.155; % eV/GPa; pressure coefficient of band gap (Booth 1982, LaBoe Vol. 22a, page 124)
% Experimental band gap as a function of temperature
EBE0=0.261; % eV
aBE0=0.026; % eV
OmegaE0=0.0189; % eV. This one is wrong in Melissa Rivero's paper.
% Full temperature dependence of the direct band gap using Bose-Einstein expression.
% This gap E0BE includes thermal expansion and electron-phonon contributions. 
E0BE=@(T)EBE0 - aBE0* (1 + 2./(exp(OmegaE0./(kB*T))-1)); % E0-exp(T) in eV
% This gap E0BEep is the pure electron-phonon contribution to the band gap.
EBE0ep=0.242; % eV
aBE0ep=0.0073; % eV
OmegaE0ep=0.0109; % eV
E0BEep=@(T)EBE0ep - aBE0ep* (1 + 2./(exp(OmegaE0ep./(kB*T))-1)); % E0-elph(T) in eV
%
% fileid=fopen('Eelph.txt','w');
E0=EBE0-aBE0; % experimental gap at 0 K
Tstart=10;
Tstep=10;
%
for i=1:1:80
    Tvar(i)=Tstart+(i-1)*Tstep;
%
    E0TE(i)=E0-3*B*dEdp*alphaint8(i)*1E-6;
%    E0TE8(i)=E0-3*B*dEdp*alphaint8(i)*1E-6;
    E0u(i)=EBE0ep-3*B*dEdp*alphaint8(i)*1E-6;
    Elph(i)=E0BE(Tvar(i))+3*B*dEdp*alphaint8(i)*1E-6;
%    
%   fprintf('T=%4.0f, alpha=%5.3f, E0(T)=%5.3f\n', ...
%   Tvar(i),alphavar(i),E0TE(i));
%    fprintf(fileID,'%5.0f, %5.3f\n',Tvar(i),Elph(i));
%    fprintf('%5.0f, %5.3f\n',Tvar(i),Elph(i));
end
%
% fclose(fileid); 
f=figure('Name','Band gap change of InSb due to thermal expansion', ...
    'Color','white');
pl=plot(Tvar,E0BE(Tvar),'--b',Tvar,E0TE,':r', ...
    Tvar,Elph,'-.g',Tvar,E0u,'-k');
%
pl(1).LineWidth=1; pl(2).LineWidth=1; pl(3).LineWidth=1; pl(4).LineWidth=2;
xlabel('Temperature T (K)',FontSize=12);
ylabel('E_0 (eV)',FontSize=12); 
ax=gca; ax.XAxis.FontSize=12; ax.YAxis.FontSize=12;
ax.XMinorTick="on"; ax.YMinorTick="on";
xlim([0 800]);
legend('experiment E_0^{exp}','thermal expansion E_0^{TE}', ...
       'electron-phonon E_0^{el-ph}','unrenormalized E_0^u', ...
       'Location','southwest','FontSize',14); 
legend('boxoff');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

