% Calculate dielectric function for interband optical transitions 
% in germanium near the direct band gap E0, following 
% Emminger et al., J. Appl. Phys. 131, 165701 (2022) and
% Menendez et al., Phys. Rev. B 101, 195204 (2020)
% Prepared by Stefan Zollner, New Mexico State University, Las Cruces, NM
% Date: 25 June 2023
% Required package: fdigamma.m 
% to calculate the digamma function with complex arguments.
% The built-in psi() function only takes real nonnegative arguments. 
%
%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 vacuum permittivity
c=299792458; % m/s speed of light in vacuum
% calculated physical constants 
h=hbar*2*pi; % eV*s
hJoule=h*e; % J*s
hbarc=hbar*c*100; % eV*cm
Rydberg=e^4*m0/(8*eps0^2*hJoule^2)/e; % eV, ground state energy for hydrogen
aBohr=(hJoule^2*eps0)/(pi*e^2*m0)*1E10; % Angstrom, Bohr radius for hydrogen
%
% amplitude prefactor for optical interband transitions
Afactor=4*sqrt(Rydberg)/3;
% optical constants n and alpha calculated from epsilon (Fox 2010, Sec. 1.3)
nc=@(z)sqrt(z); % calculates the refractive index from complex epsilon
nr=@(z)real(nc(z)); % real part of refractive index
kappa=@(z)imag(nc(z)); % extinction coefficient, imaginary part of nc
alpha=@(z,E)2*kappa(z).*E/hbarc; % absorption coefficient
%
%Material parameters for germanium (masses, matrix elements, gaps)
%Masses are at 4.2 K, from Table II in Emminger
me=0.036;
mlh=0.042;
mhh=0.33;
EP=26.0; % eV, optical matrix element at 4K, from Emminger 2022, below Eq. (22), similar to Lawaetz PRB 1971
E0=0.889; % eV, Ge direct band gap at 4K, from Emminger 2022, below Eq. (22), source: Macfarlane, Proc. Phys. Soc. 1958
% E0=0.803; % eV, Ge direct band gap at 300K, from Menendez Fig. 5b
Eu=15.6; % eV, Emminger 2022, Eq. (29), needed to calculate eps_infty
EPennT=@(T)4.146-0.05*(2/(exp(217/T)-1)+1); % eV, Penn gap, Emminger Eq. (30)
EPenn=EPennT(4); % eV, Penn gap at 4 K
eps_inftyT=@(T)1+(Eu/EPennT(T))^2; % high-frequency dielectric constant (same as static for Ge)
eps_infty=eps_inftyT(4); % Emminger 2022, Eq. (29), high-frequency dielectric constant at 4 K
%Calculated parameters
muhh=mhh*me/(mhh+me); % reduced (optical) mass, compare Emminger 2022, Table II
mulh=mlh*me/(mlh+me);
Rhh=muhh/eps_infty^2*Rydberg; % eV, exciton binding energies, compare Emminger 2022, Eq. (10)
Rlh=mulh/eps_infty^2*Rydberg; % eV
ahh=eps_infty/muhh*aBohr; % Bohr radius in Angstrom
alh=eps_infty/mulh*aBohr; % Bohr radius in Angstrom
%
Ahh=Afactor*EP*muhh^1.5; % hh to electron amplitude, compare Emminger 2022, Table II
Alh=Afactor*EP*mulh^1.5; % lh to electron amplitude, compare Emminger 2022, Table II
Gamma=0.002; % eV, broadening for transitions, compare Emminger 2022, Fig. 6
%
energy=0.80:0.001:0.98; % energy range for Emminger 2022 Fig. 1
%
% Dielectric function contribution for interband transitions (without excitons)
epshh=@(x)Ahh*(2*sqrt(E0)-sqrt(E0-x-1i*Gamma)-sqrt(x+E0+1i*Gamma))./(x+1i*Gamma).^2;
epslh=@(x)Alh*(2*sqrt(E0)-sqrt(E0-x-1i*Gamma)-sqrt(x+E0+1i*Gamma))./(x+1i*Gamma).^2;
figure('Name','Interband Transitions');
plot(energy,real(epshh(energy)+epslh(energy)),energy,imag(epshh(energy)+epslh(energy)));
title('Interband transitions at E_0 gap of Ge (no excitonic effects)');
xlabel('Energy (eV)'); ylabel('epsilon');
legend('real(epsilon)','imag(epsilon)','Location','southwest'); legend('boxoff');
%
% Dielectric function for unscreened excitonic transitions
% from C. Tanguy, Phys. Rev. Lett. 75, 4090 (1995)
% with erratum Phys. Rev. Lett. 76, 716 (1995).
xihh=@(z)sqrt(Rhh./(E0-z)); % compare Eq. (6)
xilh=@(z)sqrt(Rlh./(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)
epshhx=@(x)Ahh*sqrt(Rhh) ./(x+1i*Gamma).^2 ...
    .* ( ga(xihh(x+1i*Gamma)) + ga(xihh(-x-1i*Gamma)) -2*ga(xihh(0)) );
epslhx=@(x)Alh*sqrt(Rlh) ./(x+1i*Gamma).^2 ...
    .* ( ga(xilh(x+1i*Gamma)) + ga(xilh(-x-1i*Gamma)) -2*ga(xilh(0)) ) ;
figure('Name','Unscreened Exciton Transitions');
plot(energy,real(epshhx(energy)+epslhx(energy)),energy,imag(epshhx(energy)+epslhx(energy)));
title('Excitonic transitions at E_0 gap of Ge (unscreened)');
xlabel('Energy (eV)'); ylabel('epsilon');
legend('real(epsilon)','imag(epsilon)','Location','southwest'); legend('boxoff');
% 
% Dielectric function for screened excitonic transitions
% from C. Tanguy, Phys. Rev. B 60, 10660 (1999).
g=0.1; % select a Hulthen parameter where screening is noticeable.
xihhs=@(z)2./( sqrt((E0-z)/Rhh) + sqrt((E0-z)/Rhh+4/g) ); % compare Tanguy Eq. (7)
xilhs=@(z)2./( sqrt((E0-z)/Rlh) + sqrt((E0-z)/Rlh+4/g) ); 
gt=@(xi)-2*fdigamma(g./xi)-xi/g-2*fdigamma(1-xi)-1./xi; % compare Tanguy Eq. (8)
% The argument of epshhs and epslhs is the photon energy.
epshhs=@(x)Ahh*sqrt(Rhh)./(x+1i*Gamma).^2 ...
    .*( gt(xihhs(x+1i*Gamma)) + gt(xihhs(-x-1i*Gamma)) - 2*gt(xihhs(0)) );
epslhs=@(x)Alh*sqrt(Rlh)./(x+1i*Gamma).^2 ...
    .*( gt(xilhs(x+1i*Gamma)) + gt(xilhs(-x-1i*Gamma)) - 2*gt(xilhs(0)) );
figure('Name','Screened Exciton Transitions')
plot(energy,real(epshhs(energy)+epslhs(energy)),energy,imag(epshhs(energy)+epslhs(energy)));
title('Screened excitonic transitions at E_0 gap of Ge (g='+string(g)+')');
xlabel('Energy (eV)'); ylabel('epsilon');
legend('real(epsilon)','imag(epsilon)','Location','southwest'); legend('boxoff');
%
% Create FIGURE 1 in Emminger 2022 
g=35; % Hulthen screening parameter for undoped Ge (Menendez 2020 Fig. 5a)
Ahh=1; Rhh=0.0015; Gamma=0.001; % parameters used for Emminger Fig. 1
% These functions need to be redefined because the parameters have changed.
xihhs=@(z)2./( sqrt((E0-z)/Rhh) + sqrt((E0-z)/Rhh+4/g) ); 
xilhs=@(z)2./( sqrt((E0-z)/Rlh) + sqrt((E0-z)/Rlh+4/g) );
betas=@(x)Ahh*sqrt(Rhh) ./ (x+1i*Gamma).^2;
z1s=@(x)x+1i*Gamma;
z2s=@(x)-x-1i*Gamma;
% Unfortunately, there is an typo in the definition of the functions f1-f4.
% A factor (-2) is missing in the last term for all four functions f1-f4. 
% That error has been corrected in this program. 
% The first term describes the unbound exciton continuum response (Sommerfeld enhancement). 
f1s=@(x)-2*betas(x).* ...
    ( fdigamma(g./xihhs(z1s(x))) + fdigamma(g./xihhs(z2s(x))) - 2*fdigamma(g./xihhs(0)) );
% f2 vanishes in the absence of screening (large g)
f2s=@(x)-betas(x)/g.* (xihhs(z1s(x))+xihhs(z2s(x))-2*xihhs(0));
% The singularities of the digamma function in f3 describe the excitonic peaks.
f3s=@(x)-2*betas(x).* ...
    (fdigamma(1-xihhs(z1s(x)))+fdigamma(1-xihhs(z2s(x)))-2*fdigamma(1-xihhs(0)));
% Band-to-band transitions (no excitonic effects) are described by f4.
f4s=@(x)-betas(x).* (1./xihhs(z1s(x)) + 1./xihhs(z2s(x)) - 2./xihhs(0));
falls=@(x)f1s(x)+f2s(x)+f3s(x)+f4s(x); 
%
figure('Name','Emminger 2022 Figure 1');
subplot(2,1,1);
plot(energy,real(f1s(energy)),energy,real(f2s(energy)), ...
     energy,real(f3s(energy)),energy,real(f4s(energy)), ...
     energy,real(falls(energy)));
title('Screened excitonic transitions at E_0 gap of Ge (g='+string(g)+')');
xlabel('Energy (eV)'); ylabel('real(epsilon)');
legend('f1','f2','f3','f4','sum','Location','northwest'); legend('boxoff');
subplot(2,1,2);
plot(energy,imag(f1s(energy)),energy,imag(f2s(energy)), ...
    energy,imag(f3s(energy)),energy,imag(f4s(energy)), ...
    energy,imag(falls(energy)));
xlabel('Energy (eV)'); ylabel('imag(epsilon)');
legend('f1','f2','f3','f4','sum','Location','northwest'); legend('boxoff');
%
% Produce Fig. 5b in Menendez 2020
Gamma=0.001; % eV, broadening assumed by Menendez 2020 for Fig. 5b
energy=0.75:0.001:0.95; % energy range for Menendez 2020 Fig. 5b 
E0=0.803; % room-temperature band gap for this figure. 
% The exact results for alpha depend on the choice of epsilon_1,
% which has been assumed constant here (eps_infty).
% These functions need to be redefined because the parameters have changed.
g=100;
xihhs=@(z)2./( sqrt((E0-z)/Rhh) + sqrt((E0-z)/Rhh+4/g) ); % compare Eq. (7)
xilhs=@(z)2./( sqrt((E0-z)/Rlh) + sqrt((E0-z)/Rlh+4/g) ); 
gt=@(xi)-2*fdigamma(g./xi)-xi/g-2*fdigamma(1-xi)-1./xi; % compare Eq. (8)
epshhs1=@(x)Ahh*sqrt(Rhh)./(x+1i*Gamma).^2 ...
    .*( gt(xihhs(x+1i*Gamma)) + gt(xihhs(-x-1i*Gamma)) - 2*gt(xihhs(0)) );
epslhs1=@(x)Alh*sqrt(Rlh)./(x+1i*Gamma).^2 ...
    .*( gt(xilhs(x+1i*Gamma)) + gt(xilhs(-x-1i*Gamma)) - 2*gt(xilhs(0)) );
% These functions need to be redefined because the parameters have changed.
g=1;
xihhs=@(z)2./( sqrt((E0-z)/Rhh) + sqrt((E0-z)/Rhh+4/g) ); % compare Eq. (7)
xilhs=@(z)2./( sqrt((E0-z)/Rlh) + sqrt((E0-z)/Rlh+4/g) ); 
gt=@(xi)-2*fdigamma(g./xi)-xi/g-2*fdigamma(1-xi)-1./xi; % compare Eq. (8)
epshhs2=@(x)Ahh*sqrt(Rhh)./(x+1i*Gamma).^2 ...
    .*( gt(xihhs(x+1i*Gamma)) + gt(xihhs(-x-1i*Gamma)) - 2*gt(xihhs(0)) );
epslhs2=@(x)Alh*sqrt(Rlh)./(x+1i*Gamma).^2 ...
    .*( gt(xilhs(x+1i*Gamma)) + gt(xilhs(-x-1i*Gamma)) - 2*gt(xilhs(0)) );
% These functions need to be redefined because the parameters have changed.
g=0.3;
xihhs=@(z)2./( sqrt((E0-z)/Rhh) + sqrt((E0-z)/Rhh+4/g) ); % compare Eq. (7)
xilhs=@(z)2./( sqrt((E0-z)/Rlh) + sqrt((E0-z)/Rlh+4/g) ); 
gt=@(xi)-2*fdigamma(g./xi)-xi/g-2*fdigamma(1-xi)-1./xi; % compare Eq. (8)
epshhs3=@(x)Ahh*sqrt(Rhh)./(x+1i*Gamma).^2 ...
    .*( gt(xihhs(x+1i*Gamma)) + gt(xihhs(-x-1i*Gamma)) - 2*gt(xihhs(0)) );
epslhs3=@(x)Alh*sqrt(Rlh)./(x+1i*Gamma).^2 ...
    .*( gt(xilhs(x+1i*Gamma)) + gt(xilhs(-x-1i*Gamma)) - 2*gt(xilhs(0)) );
% These functions need to be redefined because the parameters have changed.
g=0.1;
xihhs=@(z)2./( sqrt((E0-z)/Rhh) + sqrt((E0-z)/Rhh+4/g) ); % compare Eq. (7)
xilhs=@(z)2./( sqrt((E0-z)/Rlh) + sqrt((E0-z)/Rlh+4/g) ); 
gt=@(xi)-2*fdigamma(g./xi)-xi/g-2*fdigamma(1-xi)-1./xi; % compare Eq. (8)
epshhs4=@(x)Ahh*sqrt(Rhh)./(x+1i*Gamma).^2 ...
    .*( gt(xihhs(x+1i*Gamma)) + gt(xihhs(-x-1i*Gamma)) - 2*gt(xihhs(0)) );
epslhs4=@(x)Alh*sqrt(Rlh)./(x+1i*Gamma).^2 ...
    .*( gt(xilhs(x+1i*Gamma)) + gt(xilhs(-x-1i*Gamma)) - 2*gt(xilhs(0)) );
% These functions need to be redefined because the parameters have changed.
g=0.05;
xihhs=@(z)2./( sqrt((E0-z)/Rhh) + sqrt((E0-z)/Rhh+4/g) ); % compare Eq. (7)
xilhs=@(z)2./( sqrt((E0-z)/Rlh) + sqrt((E0-z)/Rlh+4/g) ); 
gt=@(xi)-2*fdigamma(g./xi)-xi/g-2*fdigamma(1-xi)-1./xi; % compare Eq. (8)
epshhs5=@(x)Ahh*sqrt(Rhh)./(x+1i*Gamma).^2 ...
    .*( gt(xihhs(x+1i*Gamma)) + gt(xihhs(-x-1i*Gamma)) - 2*gt(xihhs(0)) );
epslhs5=@(x)Alh*sqrt(Rlh)./(x+1i*Gamma).^2 ...
    .*( gt(xilhs(x+1i*Gamma)) + gt(xilhs(-x-1i*Gamma)) - 2*gt(xilhs(0)) );
% These functions need to be redefined because the parameters have changed.
g=0.0001;
xihhs=@(z)2./( sqrt((E0-z)/Rhh) + sqrt((E0-z)/Rhh+4/g) ); % compare Eq. (7)
xilhs=@(z)2./( sqrt((E0-z)/Rlh) + sqrt((E0-z)/Rlh+4/g) ); 
gt=@(xi)-2*fdigamma(g./xi)-xi/g-2*fdigamma(1-xi)-1./xi; % compare Eq. (8)
epshhs6=@(x)Ahh*sqrt(Rhh)./(x+1i*Gamma).^2 ...
    .*( gt(xihhs(x+1i*Gamma)) + gt(xihhs(-x-1i*Gamma)) - 2*gt(xihhs(0)) );
epslhs6=@(x)Alh*sqrt(Rlh)./(x+1i*Gamma).^2 ...
    .*( gt(xilhs(x+1i*Gamma)) + gt(xilhs(-x-1i*Gamma)) - 2*gt(xilhs(0)) );
%
figure('Name','Menendez 2020 Figure 5b');
plot(energy,alpha(eps_infty+epshhs1(energy)+epslhs1(energy),energy), ...
     energy,alpha(eps_infty+epshhs2(energy)+epslhs2(energy),energy), ...
     energy,alpha(eps_infty+epshhs3(energy)+epslhs3(energy),energy), ...
     energy,alpha(eps_infty+epshhs4(energy)+epslhs4(energy),energy), ...
     energy,alpha(eps_infty+epshhs5(energy)+epslhs5(energy),energy), ...
     energy,alpha(eps_infty+epshhs6(energy)+epslhs6(energy),energy));
title('Screened excitonic absorption at E_0 gap of Ge');
xlabel('Energy (eV)'); ylabel('alpha (1/cm3)');
legend('g=100','g=1','g=0.3','g=0.1','g=0.05','g=0.0001','Location','northwest'); legend('boxoff');
%  