% Calculate interband optical transitions in alpha-tin, including band filling 
% (Burstein-Moss shift) and screened excitons. The nonparabolicity of the bands is 
% considered with adjusted effective masses. The k-dependence of the optical dipole 
% matrix element has been neglected. The formalism used is similar to Eq. (2) in 
% Carrasco, APL 113, 232104 (2018), but it uses Fermi-Dirac statistics instead of 
% Maxwell-Boltzman statistics.
%
% Prepared by Stefan Zollner, New Mexico State University, Las Cruces, NM
% Date: 22 January 2024
%
% This script will produce Figure 3 in the JVST B article by Zollner.
% Figure 2 can be produced by adjusting the effective masses in lines 509-514.
% (This may require also changing the scale of the y-axis.)
%
% 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 ...
%
% C. Emminger et al., J. Appl. Phys. 131, 165701 (2022).
% J. Menendez et al., Phys. Rev. B 101, 195204 (2020)
% and several articles by Christian Tanguy cited therein.
% R.A. Carrasco et al., APL 113, 232104 (2018)
%
% The Kramers-Kronig transformation is performed using the add-on kkbook2, see
% V. Lucarini et al, Tools for Data Analysis in Optics, Acoustics, Signal Processing;
% described in V. Lucarini et al, Kramers-Kronig Relations in Optical Materials Research,
% Springer, Heidelberg, 2005. 
% 
% Required package: fdigamma.m 
% to calculate the digamma function with complex arguments.
% The built-in psi() function only takes real nonnegative arguments. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% experimental data from Carrasco (2018)
eps2exp=[ ...
0.8	    5.8015; ...
0.79	5.7156; ...
0.78	5.6098; ...
0.77	5.5421; ...
0.76	5.4256; ...
0.75	5.3358; ...
0.74	5.2433; ...
0.73	5.1529; ...
0.72	5.0784; ...
0.71	5.0119; ...
0.7	    4.9211; ...
0.69	4.8241; ...
0.68	4.7543; ...
0.67	4.7002; ...
0.66	4.6058; ...
0.65	4.54; ...
0.64	4.4657; ...
0.63	4.402; ...
0.62	4.366; ...
0.61	4.3032; ...
0.6	    4.2453; ...
0.5998	3.9943; ...
0.59885	4.0086; ...
0.59789	3.9882; ...
0.59693	3.9922; ...
0.59598	4.0023; ...
0.59502	3.9845; ...
0.59406	3.9744; ...
0.59311	3.9576; ...
0.59215	3.9572; ...
0.59119	3.9522; ...
0.59024	3.9274; ...
0.59	4.2176; ...
0.58928	3.9188; ...
0.58832	3.9205; ...
0.58737	3.9574; ...
0.58641	3.9678; ...
0.58545	3.9412; ...
0.5845	3.9076; ...
0.58354	3.8992; ...
0.58258	3.9231; ...
0.58163	3.9109; ...
0.58067	3.8535; ...
0.58	4.1076; ...
0.57971	3.8357; ...
0.57876	3.8375; ...
0.5778	3.8379; ...
0.57684	3.8499; ...
0.57589	3.8236; ...
0.57493	3.8114; ...
0.57397	3.8444; ...
0.57302	3.8648; ...
0.57206	3.8428; ...
0.5711	3.8145; ...
0.57015	3.8226; ...
0.57	4.0522; ...
0.56919	3.8257; ...
0.56823	3.8152; ...
0.56728	3.808; ...
0.56632	3.8098; ...
0.56536	3.8056; ...
0.56441	3.7933; ...
0.56345	3.7925; ...
0.56249	3.7985; ...
0.56154	3.7825; ...
0.56058	3.7549; ...
0.56	3.9939; ...
0.55962	3.7766; ...
0.55867	3.7827; ...
0.55771	3.7704; ...
0.55675	3.7413; ...
0.5558	3.7264; ...
0.55484	3.7262; ...
0.55388	3.7314; ...
0.55293	3.751; ...
0.55197	3.7638; ...
0.55102	3.7253; ...
0.55006	3.6931; ...
0.55	3.9539; ...
0.5491	3.7082; ...
0.54814	3.709; ...
0.54719	3.6812; ...
0.54623	3.6795; ...
0.54528	3.6797; ...
0.54432	3.6811; ...
0.54336	3.7163; ...
0.54241	3.7302; ...
0.54145	3.7036; ...
0.54049	3.6802; ...
0.54	3.8866; ...
0.53953	3.663; ...
0.53858	3.6648; ...
0.53762	3.6734; ...
0.53667	3.6626; ...
0.53571	3.655; ...
0.53475	3.6643; ...
0.53379	3.6588; ...
0.53284	3.6461; ...
0.53188	3.6738; ...
0.53092	3.688; ...
0.53	3.8582; ...
0.52997	3.6575; ...
0.52901	3.636; ...
0.52806	3.6452; ...
0.5271	3.6493; ...
0.52614	3.6465; ...
0.52519	3.6476; ...
0.52423	3.6547; ...
0.52327	3.6487; ...
0.52232	3.6551; ...
0.52136	3.661; ...
0.5204	3.6428; ...
0.52	3.8828; ...
0.51945	3.6319; ...
0.51849	3.6347; ...
0.51753	3.6413; ...
0.51658	3.6346; ...
0.51562	3.64; ...
0.51466	3.6344; ...
0.51371	3.6388; ...
0.51275	3.6616; ...
0.51179	3.6683; ...
0.51084	3.673; ...
0.51	3.9336; ...
0.50988	3.663; ...
0.50892	3.6399; ...
0.50797	3.6538; ...
0.50701	3.6633; ...
0.50605	3.6684; ...
0.5051	3.6961; ...
0.50414	3.7247; ...
0.50318	3.7233; ...
0.50223	3.7257; ...
0.50127	3.7433; ...
0.50031	3.7429; ...
0.5	3.8594; ...
0.49936	3.7451; ...
0.4984	3.7605; ...
0.49744	3.7708; ...
0.49649	3.7854; ...
0.49553	3.7927; ...
0.49457	3.8157; ...
0.49362	3.8389; ...
0.49266	3.8538; ...
0.4917	3.8538; ...
0.49075	3.853; ...
0.48979	3.8794; ...
0.48883	3.9247; ...
0.48788	3.9557; ...
0.48692	3.9662; ...
0.48596	4.0054; ...
0.48501	4.042; ...
0.48405	4.0694; ...
0.48309	4.1052; ...
0.48214	4.1125; ...
0.48118	4.1268; ...
0.48022	4.1646; ...
0.47927	4.2064; ...
0.47831	4.2474; ...
0.47736	4.2974; ...
0.4764	4.3557; ...
0.47544	4.3927; ...
0.47449	4.4336; ...
0.47353	4.4832; ...
0.47257	4.5251; ...
0.47161	4.5932; ...
0.47066	4.6754; ...
0.4697	4.7295; ...
0.46874	4.798; ...
0.46779	4.8791; ...
0.46683	4.9393; ...
0.46587	4.9948; ...
0.46492	5.0776; ...
0.46396	5.1715; ...
0.46301	5.2634; ...
0.46205	5.3419; ...
0.46109	5.4286; ...
0.46013	5.5247; ...
0.45918	5.6165; ...
0.45822	5.7052; ...
0.45727	5.7892; ...
0.45631	5.8662; ...
0.45535	5.9554; ...
0.4544	6.0318; ...
0.45344	6.0996; ...
0.45248	6.1622; ...
0.45153	6.2175; ...
0.45057	6.2701; ...
0.44961	6.3123; ...
0.44866	6.3323; ...
0.4477	6.3309; ...
0.44674	6.3116; ...
0.44579	6.2982; ...
0.44483	6.2946; ...
0.44387	6.2537; ...
0.44292	6.1934; ...
0.44196	6.1041; ...
0.441	6.0029; ...
0.44005	5.8935; ...
0.43909	5.7863; ...
0.43813	5.6855; ...
0.43718	5.5857; ...
0.43622	5.5036; ...
0.43526	5.3889; ...
0.43431	5.2697; ...
0.43335	5.1743; ...
0.43239	5.0541; ...
0.43144	4.9307; ...
0.43048	4.8421; ...
0.42952	4.7507; ...
0.42857	4.6439; ...
0.42761	4.5414; ...
0.42665	4.4439; ...
0.4257	4.3625; ...
0.42474	4.2819; ...
0.42378	4.1882; ...
0.42283	4.0919; ...
0.42187	4.0251; ...
0.42091	3.9881; ...
0.41996	3.926; ...
0.419	3.8387; ...
0.41804	3.7541; ...
0.41709	3.6898; ...
0.41613	3.6452; ...
0.41517	3.593; ...
0.41422	3.5377; ...
0.41326	3.4721; ...
0.4123	3.4146; ...
0.41135	3.3767; ...
0.41039	3.3402; ...
0.40943	3.2828; ...
0.40848	3.2395; ...
0.40752	3.2343; ...
0.40657	3.209; ...
0.40561	3.1488; ...
0.40465	3.1021; ...
0.4037	3.0713; ...
0.40274	3.0481; ...
0.40178	3.0243; ...
0.40083	2.9974; ...
0.39987	2.9844; ...
0.39891	2.9653; ...
0.39796	2.9267; ...
0.397	2.8943; ...
0.39604	2.8713; ...
0.39509	2.8388; ...
0.39413	2.8191; ...
0.39317	2.8095; ...
0.39222	2.7817; ...
0.39126	2.7517; ...
0.3903	2.73; ...
0.38934	2.7011; ...
0.38839	2.674; ...
0.38743	2.6607; ...
0.38648	2.6397; ...
0.38552	2.6225; ...
0.38456	2.6215; ...
0.38361	2.6127; ...
0.38265	2.5853; ...
0.38169	2.5688; ...
0.38074	2.5668; ...
0.37978	2.541; ...
0.37882	2.5176; ...
0.37787	2.5297; ...
0.37691	2.5194; ...
0.37595	2.481; ...
0.375	2.4657; ...
0.37404	2.475; ...
0.37308	2.4689; ...
0.37213	2.4469; ...
0.37117	2.4345; ...
0.37021	2.4173; ...
0.36926	2.403; ...
0.3683	2.3992; ...
0.36734	2.3903; ...
0.36639	2.3781; ...
0.36543	2.3713; ...
0.36447	2.3585; ...
0.36352	2.3476; ...
0.36256	2.3511; ...
0.3616	2.3558; ...
0.36065	2.3365; ...
0.35969	2.3061; ...
0.35873	2.2982; ...
0.35778	2.306; ...
0.35682	2.3025; ...
0.35586	2.2898; ...
0.35491	2.2814; ...
0.35395	2.2735; ...
0.35299	2.2651; ...
0.35204	2.2539; ...
0.35108	2.2437; ...
0.35012	2.244; ...
0.34917	2.2479; ...
0.34821	2.2495; ...
0.34725	2.2444; ...
0.3463	2.2283; ...
0.34534	2.2069; ...
0.34438	2.2028; ...
0.34343	2.2088; ...
0.34247	2.1976; ...
0.34151	2.1859; ...
0.34056	2.1803; ...
0.3396	2.1798; ...
0.33864	2.1842; ...
0.33769	2.174; ...
0.33673	2.1532; ...
0.33577	2.1455; ...
0.33482	2.1455; ...
0.33386	2.141; ...
0.3329	2.1357; ...
0.33195	2.1335; ...
0.33099	2.1268; ...
0.33004	2.1217; ...
0.32908	2.1187; ...
0.32812	2.1192; ...
0.32717	2.1118; ...
0.32621	2.1012; ...
0.32525	2.1045; ...
0.3243	2.1079; ...
0.32334	2.0987; ...
0.32238	2.0852; ...
0.32142	2.0781; ...
0.32047	2.079; ...
0.31951	2.0789; ...
0.31855	2.0719; ...
0.3176	2.0674; ...
0.31664	2.0641; ...
0.31569	2.0626; ...
0.31473	2.0597; ...
0.31377	2.0489; ...
0.31282	2.0351; ...
0.31186	2.0342; ...
0.3109	2.0263; ...
0.30995	2.0189; ...
0.30899	2.0253; ...
0.30803	2.0298; ...
0.30708	2.0286; ...
0.30612	2.0184; ...
0.30516	2.0113; ...
0.30421	2.0186; ...
0.30325	2.0241; ...
0.30229	2.0173; ...
0.30134	2.0039; ...
0.30038	1.9979; ...
0.29942	1.9997; ...
0.29847	2.0027; ...
0.29751	2.0024; ...
0.29655	1.9966; ...
0.2956	1.9908; ...
0.29464	1.992; ...
0.29368	1.9986; ...
0.29273	2.0007; ...
0.29177	2.0018; ...
0.29081	1.9983; ...
0.28986	2.0034; ...
0.2889	2.0009; ...
0.28794	1.985; ...
0.28699	1.9808; ...
0.28603	1.9792; ...
0.28507	1.9735; ...
0.28412	1.9721; ...
0.28316	1.9756; ...
0.2822	1.9725; ...
0.28125	1.967; ...
0.28029	1.9712; ...
0.27933	1.9811; ...
0.27838	1.9788; ...
0.27742	1.9711; ...
0.27646	1.9725; ...
0.27551	1.9741; ...
0.27455	1.973; ...
0.27359	1.9747; ...
0.27264	1.9728; ...
0.27168	1.9725; ...
0.27072	1.9758; ...
0.26977	1.9758; ...
0.26881	1.9707; ...
0.26785	1.9642; ...
0.2669	1.9637; ...
0.26594	1.9677; ...
0.26498	1.9699; ...
0.26403	1.9681; ...
0.26307	1.9684; ...
0.26212	1.9694; ...
0.26116	1.9771; ...
0.2602	1.9847; ...
0.25925	1.9815; ...
0.25829	1.9841; ...
0.25733	1.989; ...
0.25637	1.9914; ...
0.25542	1.99; ...
0.25446	1.9798; ...
0.2535	1.9798; ...
0.25255	1.9917; ...
0.25159	1.9916; ...
0.25063	1.99; ...
0.24968	1.9934; ...
0.24872	1.9983; ...
0.24776	1.9989; ...
0.24681	1.9946; ...
0.24585	1.9905; ...
0.24489	1.995; ...
0.24394	2.0008; ...
0.24298	2.0058; ...
0.24203	2.0131; ...
0.24107	2.0143; ...
0.24011	2.0131; ...
0.23916	2.0077; ...
0.2382	2.0092; ...
0.23724	2.0268; ...
0.23629	2.0367; ...
0.23533	2.0348; ...
0.23437	2.0363; ...
0.23342	2.0339; ...
0.23246	2.0317; ...
0.2315	2.0306; ...
0.23055	2.0306; ...
0.22959	2.0332; ...
0.22863	2.0405; ...
0.22768	2.0503; ...
0.22672	2.0564; ...
0.22576	2.0493; ...
0.22481	2.0466; ...
0.22385	2.0576; ...
0.22289	2.0652; ...
0.22194	2.0594; ...
0.22098	2.0466; ...
0.22002	2.0454; ...
0.21907	2.0509; ...
0.21811	2.0452; ...
0.21715	2.0442; ...
0.2162	2.0541; ...
0.21524	2.0556; ...
0.21428	2.0474; ...
0.21333	2.0521; ...
0.21237	2.055; ...
0.21141	2.0531; ...
0.21046	2.0656; ...
0.2095	2.0592; ...
0.20854	2.0418; ...
0.20759	2.0358; ...
0.20663	2.0346; ...
0.20567	2.0298; ...
0.20472	2.0345; ...
0.20376	2.0407; ...
0.2028	2.0393; ...
0.20185	2.0326; ...
0.20089	2.0311; ...
];
%Physical constants (energy units in eV unless noted) from CODATA at NIST
%
hbar=6.582119569E-16; % eV*s reduced Planck's constant
e=1.602176634E-19; % A*s charge of the electron
m0=9.1093837015E-31; % kg free electron mass in vacuum
eps0=8.8541878128E-12; % F/m permittivity of vacuum 
c=299792458; % m/s speed of light in vacuum
kB=0.08617333262E-3; % eV/K Boltzmann constant
kBJoule=kB*e; % J Boltzmann constant in units of Joule
%
h=hbar*2*pi; % eV*s Planck's constant in eV
hJoule=h*e; % J*s Planck's constant in Joule
hbarJoule=hbar*e; % J*s reduced Planck's constant in Joule
Rydberg=e^4*m0/(8*eps0^2*hJoule^2)/e; % eV binding energy of the hydrogen atom
aBohr=(hJoule^2*eps0)/(pi*e^2*m0)*1E10; % Angstrom Bohr radius of the hydrogen atom
%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, thus factor 1E-6)
Nfactor=((2*m0*kB)/(pi*hbar*hbarJoule))^1.5*1E-6/4;
%
%Material parameters for alpha-tin (masses, matrix elements, gaps)
% m8pc=0.024; % effective mass of the \Gamma_8^+ conduction band 
m8pc=0.031; % adjusted for nonparabolicity
m8pv=0.26; % effective mass of the \Gamma_8^+ valence band (heavy hole)
% no adjustment for nonparabolicity
% m7mv=0.058; % effective mass of the \Gamma_7^- valence band
m7mv=0.1; % adjusted for nonparabolicity
% mso=0.04; % from Carrasco 2018; mso is not needed in this script
% mso=0.053; %fit to cubic characteristic equation
mL=0.19; % density of states mass for one L-valley
%
EP=24; % optical matrix element in eV
E0=0.41; % alpha-tin inverted band gap in eV (taken as a positive number)
Eind=0.005; % separation between Gamma and L point in eV
Delta0=0.8; % spin-orbit splitting
% 
mlhkp=1/((2*EP)/(-3*E0)-1);                 % effective mass from k.p theory
mekp=1/(1+EP/3*(-2/E0+1/(-E0+Delta0)));     % effective mass from k.p theory
%
eps_infty=24; % high-frequency dielectric constant (same as static)
%Calculated parameters (optical masses, exciton binding energies, exciton Bohr radii)
% optical mass for transitions from the \Gamma_7^- valence band to the \Gamma_8^+ CB
muc=(m7mv*m8pc)/(m7mv+m8pc);  
% optical mass for transitions from the \Gamma_7^- valence band to the \Gamma_8^+ VB
muv=(m7mv*m8pv)/(m8pv-m7mv);
% exciton binding energy: hole in the \Gamma_7^- VB and electron in the \Gamma_8^+ VB
Rv=muv/eps_infty^2*Rydberg; % eV
% exciton binding energy: hole in the \Gamma_7^- VB and electron in the \Gamma_8^- CB
Rc=muc/eps_infty^2*Rydberg; % eV
% Bohr radius in Angstrom: hole in the \Gamma_7^- VB and electron in the \Gamma_8^+ VB
av=eps_infty/muv*aBohr; 
% Bohr radius in Angstrom: hole in the \Gamma_7^- VB and electron in the \Gamma_8^- CB
ac=eps_infty/muc*aBohr; 
%
%Transitions into the Gamma_8^+c conduction band
Ac=Afactor*EP*muc^1.5;
%Transitions into the Gamma_8^+v valence band
Av=Afactor*EP*muv^1.5;
Gamma=0.001; % eV (broadening for transitions)
%
% Calculate Fermi energy
T=300 % K (temperature)
%
% Define Fermi-Dirac integral using polylogarithm function (Ulrich 2002)
Fermi=@(m,x)-polylog(m+1,-exp(x));
% Electron and hole densities calculated using Fermi-Dirac integrals
p=@(x)Nfactor*(m8pv*T)^1.5*Fermi(0.5,-x/(kB*T));
nL=@(x)4*Nfactor*(mL*T)^1.5*Fermi(0.5,(x-Eind)/(kB*T));
nG=@(x)Nfactor*(m8pc*T)^1.5*Fermi(0.5,x/(kB*T));
n=@(x)nL(x)+nG(x);
%
% Calculate the chemical potential mu for an intrinsic semiconductor
% where the electron density equals the hole density.
muvar=-0.1:0.01:0.1;
% This plot shows the electron and hole densities versus chemical potential
% for a given temperature T defined above. 
f=figure('Name','Electron and hole densities versus chemical potential', ...
    'Color','white');
f.Position(3:4)=[700 420];
pl=plot(muvar,n(muvar),'r-',muvar,p(muvar),'b-');
xlim([-0.1 0.1])
ax=gca; ax.FontSize=12;
xlabel('chemical potential \mu (eV)'); ylabel('carrier density (cm^{-3})'); 
legend('electron density','hole density','Location','northwest'); legend('boxoff');
pl(1).LineWidth=2; pl(2).LineWidth=2;
text(0.05,0.7,'T='+string(T)+' K','FontSize',14,'Units','normalized');
%
diff=@(x)n(x)-p(x);
muroot=fzero(diff,0);
mu=muroot;
mumeV=mu*1000
text(0.05,0.6,'\mu='+string(round(mumeV,3,"significant"))+' meV','FontSize',14,'Units','normalized'); 
%
% Carrier concentrations n and p in the valence and conduction bands (Gamma and L)
% for the chemical potential calculated above (and temperature T as defined).
%
proot=p(muroot)
text(0.05,0.5,'n=p='+string(round(proot,3,"significant"))+' cm^{-3}','FontSize',14,'Units','normalized'); 
nGroot=nG(muroot);
nLroot=nL(muroot);
close(f);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Define dielectric function contribution for transitions from Gamma_7^-v valence band.
% This is for band-to-band transitions without consideration of excitonic effects.
% The Pauli blocking term with Fermi function is not included here. 
% See Equations (6) and (14) in alpha-tin manuscript. 
%
epsc=@(x)Ac*(2*sqrt(E0)-sqrt(E0-x-1i*Gamma)-sqrt(x+E0+1i*Gamma))./(x+1i*Gamma).^2;
epsv=@(x)Av*(2*sqrt(E0)-sqrt(E0-x-1i*Gamma)-sqrt(x+E0+1i*Gamma))./(x+1i*Gamma).^2;
% Define Fermi-Dirac distribution function
f=@(x)1./(exp((x-mu)/(kB*T))+1);
% Occupation factors in Gamma8pc and Gamma8pv bands as a function of photon energy
% See Equations (10) and (17) in alpha-tin manuscript. 
f8pc=@(x)f((m7mv/(m7mv+m8pc))*(x-E0));
f8pv=@(x)f(-(m7mv/(m8pv-m7mv))*(x-E0));
%
% Calculate epsilon_2 with Pauli blocking included.
% See Equations (6) and (7) in alpha-tin manuscript.
epscFD2=@(x)imag(epsc(x)).*(1-f8pc(x));
epsvFD2=@(x)imag(epsv(x)).*(1-f8pv(x));
%
% Kramers-Kronig tests suggest that an energy range from 0 to 20 eV with 2 meV steps
% is required for an accurate KK transform to obtain the real part. 
energy=0.0:0.002:20;
% Plot Fermi occupation factors for the transition
% plot(energy,f(energy));
% plot(energy,f8pc(energy));
% plot(energy,f8pv(energy));
% xlim([0.2 0.8]);
%
% Kramers-Kronig tests (without occupation factors)
% epsc(energy);
% epsv(energy);
% epsC2=imag(epsc(energy));
% epsV2=imag(epsv(energy));
% epsC1=real(epsc(energy));
% epsV1=real(epsv(energy)); 
% epsC1=kkrebook2(energy,epsC2,0);
%
% Calculate dielectric function with occupation factors
epsC2=epscFD2(energy); % imaginary part calculated analytically
epsV2=epsvFD2(energy);
epsC1=kkrebook2(energy,epsC2,0); % real part obtained from KK transform
epsV1=kkrebook2(energy,epsV2,0);
%
% Plot the various contributions to interband transitions in inverted semiconductors.
% Compare Fig. 3 in Carrasco (2018). 
f=figure( ...
    'Name', 'Interband transitions in inverted-gap alpha-tin', ...
    'Color','white',Resize='off',Units='inches',Position=[0 0 3.378 5]);
tiledlayout("vertical","TileSpacing","none");
nexttile;
pl=plot(energy,epsV2,'b:',energy,epsC2,'g--',energy,epsC2+epsV2,'r-');
% title('Interband transitions in inverted-gap alpha-tin (no excitonic effects)');
ax=gca; ax.FontSize=10;
% xlabel('Energy (eV)'); 
ylabel('\epsilon_2'); 
xlim([0.2 0.8]);
ax.XTick=0.2:0.1:0.8;
xticklabels({}); 
% ylim([-0.5 3]); for parabolic bands
% yticklabels({'',0,0.5,1,1.5,2,2.5});
ylim([-1 8]); % for nonparabolic bands
pl(1).LineWidth=1; pl(2).LineWidth=1; pl(3).LineWidth=1; 
% legend('\Gamma_8^{+v}','\Gamma_8^{+v}','\Gamma_8^{+c}','\Gamma_8^{+c}','total','total');
text(0.05,0.9,'(a)','FontSize',10,'Units','normalized'); 
text(0.05,0.5,'T='+string(T)+' K','FontSize',10,'Units','normalized');
% add experimental data from Carrasco 2018 to the plot.
hold on;
% pl=plot(eps2exp(1:448,1),eps2exp(1:448,2),'+');
pl=plot(eps2exp(1:448,1),eps2exp(1:448,2)-2,'+'); % subtract background of 2
hold off; 
lgd=legend('\Gamma_8^{+v}','\Gamma_8^{+c}','total','exp',"Location","northeast");
lgd.FontSize=8;
legend('boxoff'); 
% text(0.05,0.8,'E_{ind}='+string(Eind*1000)+' meV','FontSize',14,'Units','normalized');
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate dielectric function for full excitonic effects without screening,
% see C. Tanguy, Phys. Rev. Lett. 75, 4090 (1995)
% with erratum Phys. Rev. Lett. 76, 716 (1995).
xiv=@(z)sqrt(Rv./(E0-z)); % compare Eq. (6) in Tanguy
xic=@(z)sqrt(Rc./(E0-z));
ga=@(z)2*log(z)-2*pi*cot(pi*z)-2*fdigamma(z)-1./z; % compare Eq. (12) of Erratum
% compare Tanguy 1995 Eq. (11)
epsvx=@(x)Av*sqrt(Rv) ./(x+1i*Gamma).^2 ...
    .* ( ga(xiv(x+1i*Gamma)) + ga(xiv(-x-1i*Gamma)) -2*ga(xiv(0)) );
epscx=@(x)Ac*sqrt(Rc) ./(x+1i*Gamma).^2 ...
    .* ( ga(xic(x+1i*Gamma)) + ga(xic(-x-1i*Gamma)) -2*ga(xic(0)) ) ;
% figure('Name','Unscreened Exciton Transitions, without occupation factors');
% plot(energy,real(epsvx(energy)+epscx(energy)),energy,imag(epsvx(energy)+epscx(energy)));
% title('Excitonic transitions at E0 gap of a-tin (unscreened, no occupation factors)');
%
% Now include occupation factors
epscFD2x=@(x)imag(epscx(x)).*(1-f8pc(x));
epsvFD2x=@(x)imag(epsvx(x)).*(1-f8pv(x));
epsC2x=epscFD2x(energy);
epsV2x=epsvFD2x(energy);
% Do the Kramers-Kronig transform
epsC1x=kkrebook2(energy,epsC2x,0);
epsV1x=kkrebook2(energy,epsV2x,0);
% figure('Name','Unscreened Exciton Transitions, with occupation factors');
% pl=plot(energy,epsV1x,'b:',energy,epsV2x,'b:', ...
%     energy,epsC1x,'g--',energy,epsC2x,'g--', ...
%     energy,epsC1x+epsV1x,'r-',energy,epsC2x+epsV2x,'r-');
% pl(1).LineWidth=1; pl(2).LineWidth=3; pl(3).LineWidth=1; pl(4).LineWidth=3;
% pl(5).LineWidth=1; pl(6).LineWidth=3;
% title('Excitonic transitions at E0 gap of a-tin (unscreened, with occupation factors)');
% xlabel('Energy (eV)'); ylabel('epsilon'); xlim([0.2 0.8]);
% legend('real(epsilon)','imag(epsilon)','Location','southwest'); legend('boxoff');
% ylim([-1 3]);
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Screening parameters: 
% Debye wave vector kD (in inverse cm), kDA (in inverse Angstrom)
% Debye screen length lambdaDA, Hulthen parameters gv and gc for G8+V and G8+C excitons.
% Only consider screening by holes (density proot), 
% because most electrons are at the L point (at room temperature). 
%
kD=sqrt(DebyeFactor*proot/(eps_infty*T)); % inverse cm
kDA=kD*1E-8; % inverse Angstrom
lambdaDA=1/kDA; % Angstrom
gv=12/(pi^2*av*kDA)
gc=12/(pi^2*ac*kDA)
%
%  screened exciton absorption from Tanguy 1999
%
xicg=@(z) 2 ./ ( sqrt((E0-z)/Rc) + sqrt((E0-z)/Rc+4/gc) );
xivg=@(z) 2 ./ ( sqrt((E0-z)/Rv) + sqrt((E0-z)/Rv+4/gv) );
gtc=@(z) -2*fdigamma(gc./z)-z/gc-2*fdigamma(1-z)-1./z;
gtv=@(z) -2*fdigamma(gv./z)-z/gv-2*fdigamma(1-z)-1./z;
%
epsvg=@(x)Av*sqrt(Rv) ./(x+1i*Gamma).^2 ...
    .* ( gtv(xivg(x+1i*Gamma)) + gtv(xivg(-x-1i*Gamma)) -2*gtv(xivg(0)) );
epscg=@(x)Ac*sqrt(Rc) ./(x+1i*Gamma).^2 ...
    .* ( gtc(xicg(x+1i*Gamma)) + gtc(xicg(-x-1i*Gamma)) -2*gtc(xicg(0)) ) ;
%
epscFD2g=@(x)imag(epscg(x)).*(1-f8pc(x));
epsvFD2g=@(x)imag(epsvg(x)).*(1-f8pv(x));
epsC2g=epscFD2g(energy);
epsV2g=epsvFD2g(energy);
% Do the Kramers-Kronig transform 
epsC1g=kkrebook2(energy,epsC2g,0);
epsV1g=kkrebook2(energy,epsV2g,0);
%
% f=figure('Name','Screened Exciton Transitions, with occupation factors', ...
%         'Color','white',Resize='off',Units='inches',Position=[0 0 3.378 3]);
nexttile;
% pl=plot(energy,epsV1g,'b:',energy,epsV2g,'b:', ...
%     energy,epsC1g,'g--',energy,epsC2g,'g--', ...
%     energy,epsC1g+epsV1g,'r-',energy,epsC2g+epsV2g,'r-');
pl=plot(energy,epsV2g,'b:',energy,epsC2g,'g--',energy,epsC2g+epsV2g,'r-');
% title('Excitonic transitions at E0 gap of a-tin (screened, with occupation factors)');
ax=gca; ax.FontSize=10;
xlabel('Energy (eV)'); 
lgd=legend('\Gamma_8^{+v}','\Gamma_8^{+c}','total');
lgd.FontSize=8;
ylabel('\epsilon_2');
% ylim([-0.5 3]);  % for parabolic bands
% yticklabels({'',0,0.5,1,1.5,2,2.5});
ylim([-1 8]); % for nonparabolic bands
pl(1).LineWidth=1; pl(2).LineWidth=1; pl(3).LineWidth=1;
% pl(4).LineWidth=2; pl(5).LineWidth=1; pl(6).LineWidth=2;
% legend('\Gamma_8^{+v}','\Gamma_8^{+v}','\Gamma_8^{+c}','\Gamma_8^{+c}','total','total');
text(0.05,0.5,'T='+string(T)+' K','FontSize',10,'Units','normalized');
text(0.05,0.9,'(b)','FontSize',10,'Units','normalized'); 
xlim([0.2 0.8]); ax.XTick=0.2:0.1:0.8;
% text(0.05,0.8,'E_{ind}='+string(Eind*1000)+' meV','FontSize',14,'Units','normalized');
% return;
% add experimental data from Carrasco 2018 to the plot.
hold on;
% pl=plot(eps2exp(1:448,1),eps2exp(1:448,2),'+');
pl=plot(eps2exp(1:448,1),eps2exp(1:448,2)-2,'+'); % subtract background of 2
hold off; 
legend('\Gamma_8^{+v}','\Gamma_8^{+c}','total','exp',"Location","northeast");
legend('boxoff'); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

