% Band structure of alpha-tin following Kane (1957) including non-parabolicity 
% corrections, calculated from Vieta's solution to the cubic equation.
%
% Prepared by Stefan Zollner, New Mexico State University, Las Cruces, NM
% Date: 22 January 2024
%
% This script will produce Figure 1 in the JVST B article by Zollner.
%
% References:
%
% "Excitonic Effects in the Optical Absorption of Gapless Semiconductor alpha-Tin
% Near the Direct Band Gap", accepted by Journal of Vacuum Science and Technology B  
% on 01/05/2024, published on ...
% 
% Evan O. Kane, "Band structure of indium antimonide", 
% J. Phys. Chem. Solids 1, p. 249 (1957).
%
%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
c=299792458; % m/s speed of light in vacuum
%
kBJoule=kB*e;
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
%
% Material parameters for alpha-tin (compare Ioffe Institute web site)
% a=???; % lattice constant of alpha-tin in Angstrom, probably 300 K
E0=-0.413; % eV, from Lawaetz 1971, at 0 K
Delta0=0.8; % eV
EP=23.8; % eV, momentum matrix energy, from Lawaetz 1971, at 0 K
% EP=0.44*Hartree*2; % convert 2P^2 from atomic units to eV, from Kane 1957
% me=1/(1+EP/3*(2/E0+1/(E0+Delta0))); % should be about -0.058 (Lawaetz)
mhh=0.26; % from Carrasco 2018
me=-0.058; % from Carrasco 2018
% me=-0.1; %fit to cubic characteristic equation
% mlh=me; % (within this 4x4 k.p model)
mlh=-0.024; % from Carrasco 2018
% mlh=-0.031; %fit to cubic characteristic equation
mso=0.04; % from Carrasco 2018
% mso=0.053; %fit to cubic characteristic equation
%
% Dispersion of parabolic bands calculated with masses specified above. 
Ehh=@(k)hbar2m0*k.^2/2; % heavy hole from Kane's 8x8 Hamiltonian
Ecp=@(k)E0+hbar2m0*k.^2/(2*me);
Elh=@(k)-hbar2m0*k.^2/(2*mlh);
Ehhp=@(k)-hbar2m0*k.^2/(2*mhh);
Esop=@(k)-Delta0-hbar2m0*k.^2/(2*mso);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Vieta's solution to the cubic characteristic equation
k=0:0.001:0.1; % wave vector k for band structure plot
kau=k*aBohr; % wave vector k in atomic units
%
Asqu=@(k)hbar2m0*k.^2/2*EP/3;
a=1;
b=Delta0-E0;
c=@(k)-E0*Delta0-3*Asqu(k);
d=@(k)-2*Asqu(k)*Delta0;
% use a=1
p=@(k)c(k)-b^2/3;
q=@(k)2*b^3/27-b*c(k)/3+d(k);
%
t0=@(k)2*sqrt(-p(k)/3).* ...
    cos( acos( 3*q(k)/2./p(k) .*sqrt(-3./p(k))) /3 );
t1=@(k)2*sqrt(-p(k)/3).* ...
    cos( acos( 3*q(k)/2./p(k) .*sqrt(-3./p(k))) /3 -2*pi/3);
t2=@(k)2*sqrt(-p(k)/3).* ...
    cos( acos( 3*q(k)/2./p(k) .*sqrt(-3./p(k))) /3 -4*pi/3);
% 
Ene0=@(k)t0(k) -(Delta0-E0)/3 +hbar2m0*k.^2/2;
Ene1=@(k)t1(k) -(Delta0-E0)/3 +hbar2m0*k.^2/2;
Ene2=@(k)t2(k) -(Delta0-E0)/3 +hbar2m0*k.^2/2;
%
% Draw a figure with the band structure, compare 8x8 k.p and parabolic bands.
%
fig=figure(Name='Bands in alpha-Sn from 8x8 k.p model with Vietas equation', ...
    Color='white',Resize='off',Units='inches',Position=[12 3 3.378 4]);
%
pl=plot(kau,Ene0(k),'r',kau,Ene1(k),'g',kau,Ene2(k),'k', ... % kau,Ehh(k),'b', ...
        kau,Ecp(k),'g--',kau,Elh(k),'r--',kau,Esop(k),'k--',kau,Ehhp(k),'b--');
pl(1).LineWidth=1; pl(2).LineWidth=1; pl(3).LineWidth=1; % pl(4).LineWidth=1; 
pl(4).LineWidth=1; pl(5).LineWidth=1; pl(6).LineWidth=1; pl(7).LineWidth=1;
ax=gca; ax.FontSize=9; %
xlim([0,0.02]);
ax.XTick=0.0:0.005:0.02; ax.XMinorTick="on";
ylim([-0.9,0.2]);
ax.YTick=-0.9:0.1:0.2; ax.YMinorTick="on";
EFermi=-0.012;  % Fermi energy for a-tin at 300 K
yl=yline(EFermi,'k-','E_F','Linewidth',1); % draw a line at the Fermi level
yl.LabelHorizontalAlignment='left';
yl.LabelVerticalAlignment='bottom';
yl.FontSize=10;
annotation('arrow',[0.5 0.5],[0.46,0.77],'Color','b','LineWidth',1, ...
    HeadLength=5,HeadWidth=5);
annotation('arrow',[0.8 0.8],[0.42,0.90],'Color','r','LineWidth',1, ...
    HeadLength=5,HeadWidth=5);
annotation('arrow',[0.6 0.6],[0.78,0.84],'Color','k','LineWidth',1,'LineStyle',':', ...
    HeadLength=5,HeadWidth=5);
annotation('arrow',[0.55 0.55],[0.16,0.83],'Color','k','LineWidth',1,'LineStyle',':', ...
    HeadLength=5,HeadWidth=5);
annotation('arrow',[0.575 0.575],[0.15,0.77],'Color','k','LineWidth',1,'LineStyle',':', ...
    HeadLength=5,HeadWidth=5);
%
xlabel('wave vector k (atomic units)'); ylabel('Energy (eV)');
text(0.0005,-0.7,'\Gamma_7^+ (so)','FontSize',9,'Units','data','Color','k');
text(0.005,-0.50,'\Gamma_7^- (e)','FontSize',9,'Units','data','Color','g');
text(0.005,0.12,'\Gamma_8^+ (lh)','FontSize',9,'Units','data','Color','r');
text(0.012,-0.10,'\Gamma_8^+ (hh)','FontSize',9,'Units','data','Color','b');
text(0.014,-0.7,'k_F','FontSize',9,'Units','data','Color','k');
%
% end
%

