In [1]:
# Hellwarth 1999 PRB - Part IV; T-dep of the Feynman variation parameter
# A Friday afternoon of hacking to try and implement the T-dep electron-phonon coupling from the above PRB
# Which was unusually successful! And more or less reproduced Table III
In [2]:
# Just in case anyone is following this from the far future; we are using:
versioninfo()
In [3]:
# one-dimensional numerical integration in Julia using adaptive Gauss-Kronrod quadrature
using QuadGK
# Equation numbers follow above Hellwarth 1999 PRB
# 62b
A(v,w,β)=3/β*( log(v/w) - 1/2*log(2*π*β) - log(sinh(v*β/2)/sinh(w*β/2)))
# 62d
Y(x,v,β)=1/(1-exp(-v*β))*(1+exp(-v*β)-exp(-v*x)-exp(v*(x-β)))
# 62c integrand
f(x,v,w,β)=(exp(β-x)+exp(x))/(w^2*x*(1-x/β)+Y(x,v,β)*(v^2-w^2)/v)^(1/2)
# 62c
B(v,w,β,α) = α*v/(sqrt(π)*(exp(β)-1)) * quadgk(x->f(x,v,w,β),0,β/2)[1]
#62e
C(v,w,β)=3/4*(v^2-w^2)/v * (coth(v*β/2)-2/(v*β))
F(v,w,β,α)=-(A(v,w,β)+B(v,w,β,α)+C(v,w,β)) #(62a)
# Can now evaluate, e.g.
# F(v,w,β,α)=F(7.2,6.5,1.0,1.0)
# BUT - this is just the objective function! Not the optimised parameters.
# Also there's a scary numeric integration (quadgk) buried within...
Out[3]:
In [4]:
# Print out F(alpha,beta) for a specific v,w; as a test
@printf("\t\t")
for α in 1:5
@printf("α=%d\t\t",α)
end
@printf("\n")
for β in 1:0.25:3.0
v=w=4
print("β: $β \t||")
for α in 1:5
@printf("%f\t",F(v,w,β,α))
end
println()
end
In [5]:
# OK - very primitive!
# But these are 1D traces along the solution for Alpha=Beta=1 in Helwarth PRB TABLE III,
# this was used to correct a transcription error in the above typed-in equations
# It was also good to see what F(v,w) looked like as a function of v and w near an optimal solution
v=7.20
w=6.5
α=1.0
β=1.0
for v=6:0.1:8
@printf("%f %f\n",v,F(v,w,β,α))
end
@printf("\n")
v=7.20
for w=6:0.1:7
@printf("%f %f\n",w,F(v,w,β,α))
end
In [6]:
# Angle for the ringside seats, when the fall, don't blame me, Bring on the Major Leagues
using Optim
# Julia package stuffed full of magic, does auto-differentation & etc. etc.
In [ ]:
In [7]:
Fopt(x) = F(x[1],x[2],1,1)
Fopt([7.2,6.5])
# OK! It looks like I can bury the alpha, beta parameters (which we don't optimise), by wrapping our function in a function definition.
Out[7]:
In [8]:
initial=[7.2,6.5]
optimize(Fopt, initial, LBFGS())
Out[8]:
In [9]:
optimize(Fopt, initial, BFGS(), Optim.Options(autodiff=true))
Out[9]:
In [10]:
optimize(Fopt, initial, BFGS(), Optim.Options(autodiff=true))
Out[10]:
In [11]:
# Right! The above looks like this might actually just work...
@printf("\t\t")
for α in 1:5
@printf("α=%d\t\t",α)
end
@printf("\n")
for β in 1:0.25:3.0
print("β: $β \t||")
for α in 1:5
myf(x) = F(x[1],x[2],β,α)
solution=Optim.minimizer(optimize(myf, initial,ConjugateGradient(), Optim.Options(autodiff=true)))
#print(solution,"\t")
@printf("%.2f %.2f\t",solution[1],solution[2])
end
println()
end
In [12]:
# After a bit of fiddling, I figured out how to add bounds, to stop that 'DomainError',
# which occurs where the you are evaluating log(-ve Real), i.e. w<0.0 or v<0.0
initial=[7.2,6.5]
lower=[0.0,0.0]
upper=[10.0,10.0]
@printf("\t\t")
for α in 1:5
@printf("α=%d\t\t",α)
end
@printf("\n")
for β in 1:0.25:3.0
print("β: $β \t||")
for α in 1:5
myf(x) = F(x[1],x[2],β,α)
solution=optimize(DifferentiableFunction(myf), initial, lower, upper, Fminbox(); optimizer = ConjugateGradient, optimizer_o=Optim.Options(autodiff=true))
minimum=Optim.minimizer(solution)
v=minimum[1]
w=minimum[2]
#print(solution,"\t")
@printf("%.2f %.2f\t",v,w)
end
println()
end
In [13]:
# So that looks really good! I was super stoked to see how close these values are to TABLE III in Hellwarth
# However, the solutions all start on (7.20,6.50) so that top-left data point is cheating, whereas the
# others have some disagreement / noise associated with them
# I was wondering whether it might be a function of the optimiser, so thought I'd try them all
initial=[7.1,6.5]
# Main use of these bounds is stopping v or w going negative, at which you get a NaN error as you are evaluating log(-ve Real)
lower=[1.0,1.0]
upper=[10.0,10.0]
for optimizer in [BFGS, LBFGS, ConjugateGradient] # Newton, GradientDescent, NelderMead - steps outside box & log(-ve)->NaN error
@printf("\n\t\t##### NOW TRIALING: %s #####\n\n",optimizer)
@printf("\t\t")
for α in 1:5
@printf("α=%d\t\t",α)
end
@printf("\n")
for β in 1:0.25:3.0
print("β: $β \t||")
for α in 1:5
myf(x) = F(x[1],x[2],β,α)
res=optimize(DifferentiableFunction(myf), initial, lower, upper, Fminbox(); optimizer = optimizer, optimizer_o=Optim.Options(autodiff=true))
minimum=Optim.minimizer(res)
#show(Optim.converged(res)) # All came out as 'true'
#print(solution,"\t")
@printf("%.2f %.2f\t",minimum[1],minimum[2])
end
println()
end
end
In [14]:
# So why the (slight) disagreement? I don't really know.
# It may well be that Julia has a much better control of errors,
# certainly the Optim package (optimisation) defaults seem to be at the Machine-precision end of the world,
# and a lot of intermediates will be Float64 or larger.
# The objective function seems quite flat in the vicinity of the optimum, so this would also explain some of the variance.
In [15]:
# OK, now Ref 9; Schultz PR Volume 116, Number 3, Nov 1959
const MassElectron = 9.10938188e-31; # kg
v=7.15
w=6.51
m=1.0 #MassElectron
# Eqn (2.4)
μ=m*(v^2-w^2)/v^2
println("Reduced mass: $μ")
rf=(3/2*μ*v)^(1/2)
println("Feynman Polaron Radius: $rf")
# Units?!? Bohr?
In [16]:
using Optim
const hbar = const ħ = 1.05457162825e-34; # kg m2 / s
const eV = const q = const ElectronVolt = 1.602176487e-19; # kg m2 / s2
const me=MassElectron = 9.10938188e-31; # kg
const Boltzmann = const kB = 1.3806504e-23; # kg m2 / K s2
# "Band Mass" in Deveresse, actually m_b
mb=0.15*me #*1.4 # 0.15 is band reduced mass; 1.4 is polaron renormalised mass
freq=2.25E12
ω = (2*pi)*freq
α=2.395939683378253
@printf("α=%f\n",α)
initial=[7.1,6.5]
#initial=[9.0,9.0]
# Main use of these bounds is stopping v or w going negative, at which you get a NaN error as you are evaluating log(-ve Real)
lower=[0.1,0.1]
upper=[100.0,100.0]
# Empty arrays for storing data
Ts=[]
Kμs=[]
Hμs=[]
FHIPμs=[]
ks=[]
Ms=[]
As=[]
Bs=[]
Cs=[]
Fs=[]
# I think this β term subsumes the energy of the phonon; i.e. kbT c.f. ħω
for T in 10:10:300
β=1/(kB*T)
βred=ħ*ω*β
@printf("T: %f β: %.2g βred: %.2g\t",T,β,βred)
myf(x) = F(x[1],x[2],βred,α)
res=optimize(DifferentiableFunction(myf), initial, lower, upper, Fminbox(); optimizer = BFGS, optimizer_o=Optim.Options(autodiff=true))
minimum=Optim.minimizer(res)
#show(Optim.converged(res)) # All came out as 'true'
v=minimum[1]
w=minimum[2]
@printf(" v= %.2f w= %.2f\t",v,w)
# From 1962 Feynman, definition of v and w in terms of the coupled Mass and spring-constant
# See Page 1007, just after equation (18)
# Units of M appear to be 'electron masses'
# Unsure of units for k, spring coupling constant
k=(v^2-w^2)
M=(v^2-w^2)/w^2
@printf(" M=%f k=%f\t",M,k)
append!(ks,k)
append!(Ms,M)
# F(v,w,β,α)=-(A(v,w,β)+B(v,w,β,α)+C(v,w,β)) #(62a) - Hellwarth 1999
@printf("\n Polaron Es: A= %f B= %f C= %f F= %f",A(v,w,βred),B(v,w,βred,α),C(v,w,βred),F(v,w,βred,α))
append!(As,A(v,w,βred))
append!(Bs,B(v,w,βred,α))
append!(Cs,C(v,w,βred))
append!(Fs,F(v,w,βred,α))
#[1.60] in Devereese, page 36; 6th Edition of Frohlich polaron notes
# I believe here β is in SI (expanded) units
μ=(w/v)^3 * (3*q)/(4*mb*ħ*ω^2*α*β) * exp(ħ*ω*β) * exp((v^2-w^2)/(w^2*v))
@printf("\n\tμ(FHIP)= %f m^2/Vs \t= %.2f cm^2/Vs",μ,μ*100^2)
append!(Ts,T)
append!(FHIPμs,μ*100^2)
#[1.61] - Kadanoff's Boltzmann eqn derived mob
μ=(w/v)^3 * (q)/(2*mb*ω*α) * exp(ħ*ω*β) * exp((v^2-w^2)/(w^2*v))
@printf("\n\tμ(Kadanoff)= %f m^2/Vs \t= %.2f cm^2/Vs",μ,μ*100^2)
append!(Kμs,μ*100^2)
#Hellwarth1999 Eqn (2) and (1) - Unclear where from (just references FHIP), contains factors of βred
R=(v^2-w^2)/(w^2*v) # inline, page 300 just after Eqn (2)
#b=R*βred/sinh(b*βred*v/2) # this self-references b! what the hell?
# Oops - it's a typo from Feynman1962!
b=R*βred/sinh(βred*v/2) # Feynman1962 version; page 1010, Eqn (47b)
#b=0 # Hellwarth1999/Baggio1997 "Setting b=0 makes less than 0.1% error"
a=sqrt( (βred/2)^2 + R*βred*coth(βred*v/2))
k(u,a,b,v) = (u^2+a^2-b*cos(v*u))^(-3/2)*cos(u) # integrand in (2)
K=quadgk(u->k(u,a,b,v),0,Inf)[1] # numerical quadrature integration of (2)
#Right-hand-side of Eqn 1 in Hellwarth 1999 // Eqn (4) in Baggio1997
RHS= α/(3*sqrt(π)) * βred^(5/2) / sinh(βred/2) * (v^3/w^3) * K
μ=RHS^-1 * (q)/(ω*mb)
@printf("\n\tμ(Hellwarth1999)= %f m^2/Vs \t= %.2f cm^2/Vs",μ,μ*100^2)
append!(Hμs,μ*100^2)
#Hellwarth1999 b=0
R=(v^2-w^2)/(w^2*v) # inline, page 300 just after Eqn (2)
b=0 # Hellwarth1999/Baggio1997 "Setting b=0 makes less than 0.1% error"
a=sqrt( (βred/2)^2 + R*βred*coth(βred*v/2))
k(u,a,b,v) = (u^2+a^2-b*cos(v*u))^(-3/2)*cos(u) # integrand in (2)
K=quadgk(u->k(u,a,b,v),0,Inf)[1] # numerical quadrature integration of (2)
#Right-hand-side of Eqn 1 in Hellwarth 1999 // Eqn (4) in Baggio1997
RHS= α/(3*sqrt(π)) * βred^(5/2) / sinh(βred/2) * (v^3/w^3) * K
μ=RHS^-1 * (q)/(ω*mb)
@printf("\n\tμ(Hellwarth1999,b=0)= %f m^2/Vs \t= %.2f cm^2/Vs",μ,μ*100^2)
#append!(Hμs,μ*100^2)
@printf("\n\n")
# Recycle previous results as next guess
initial=[v,w] # Caution! Might cause weird sticking in local minima
end
In [17]:
using Plots
default(size=(600,400))
In [18]:
plot(Ts,Ms,label="Mass",marker=2,xlab="Temperature (K)",ylab="Mass of Phonon cloud (electron masses)",ylim=(0,1.2))
#plot!([0],[0])
#plot!(Ts,ks,label="Spring Consts",marker=2)
Out[18]:
In [19]:
plot(Ts,ks,label="Spring Consts",marker=2, xlab="Temperature (K)",ylab="Some internal spring const",)
plot!([0],[0])
Out[19]:
In [20]:
plot(Ts,As,label="A",marker=2, xlab="Temperature (K)",ylab="Energy ?")
plot!(Ts,Bs,label="B",marker=2)
plot!(Ts,Cs,label="C",marker=2)
#plot!(Ts,Fs,label="F",marker=2)
plot!([0],[0])
Out[20]:
In [21]:
using Plots
default(size=(600,400))
plot(Ts,Kμs,label="Kadanoff Polaron mobility",marker=2,xlab="Temperature (K)",ylab="Mobility (cm^2/Vs)",ylims=(0,1000))
plot!(Ts,FHIPμs,label="FHIP",marker=2)
plot!(Ts,Hμs,label="Hellwarth1999",marker=2)
#plot!([0],[0])
Out[21]:
In [22]:
# Milot/Herz 2015 Time-Resolved-Microwave-Conductivity mobilities
# Data from table in SI of: DOI: 10.1002/adfm.201502340
# Absolute values possibly dodge due to unknown yield of charge carriers; but hopefully trend A.OK!
Milot= [
8 184
40 321
80 143
120 62
140 40
160 52
180 44
205 41
230 39
265 26
295 35
310 24
320 24
330 19
340 16
355 15
]
# IV estimated mobilities (?!) from large single crystals, assumed ambient T
# Nature Communications 6, Article number: 7586 (2015)
# doi:10.1038/ncomms8586
Saidaminov =
[ 300 67.2 ]
#Semonin2016,
# doi = {10.1021/acs.jpclett.6b01308},
Semonin =
[ 300 115 ] # +- 15 cm^2/Vs, holes+electrons
Out[22]:
In [23]:
#default(size=(1200,800))
plot(Milot[:,1],Milot[:,2],label="Milot T-dep TRMC Polycrystal",
xlab="Temperature (K)",ylab="Mobility (cm^2/Vs)",marker=2, ylims=(0,400) )
plot!(Saidaminov[:,1],Saidaminov[:,2],label="Saidaminov JV Single Crystal", marker=6)
plot!(Semonin[:,1],Semonin[:,2],label="Semonin Single Crystal TRMC", marker=6)
plot!(Ts,Kμs,label="Kadanoff Polaron mobility",marker=2)
plot!(Ts,Hμs,label="Hellwarth1999 Polaron mobility",marker=2)
#plot!([0],[0])
Out[23]:
In [24]:
plot!(yscale=(:log))
In [25]:
default(size=(1200,800))
In [26]:
# OK, now Ref 9 in Hellwarth 1999 -
# Schultz Physical Review Volume 116, Number 3, Nov 1959
# Slow Electrons in Polar Crystals: Self-Energy, Mass, and Mobility
const MassElectron = 9.10938188e-31; # kg
# 300 K fit from above
v=19.83
w= 16.93
# Guess at What Hellwarth 1999 used:
# "Using these values we estimate a room-temperature polaron radius of approx 0.6 nm for BiSiO"
# So they're using an omega=500 cm^-1 = 62 meV, so βred=2.5
# Not sure what alpha they are using; (43) gives it in terms of phonon energies + reduced masses
# Best guess for now - alpha=3, βred=2.5, so from Table III
v=4.10
w=2.55
m=1.0 #MassElectron
# Eqn (2.4), Shultz1959
μ=m*(v^2-w^2)/v^2
println("Reduced mass: $μ")
rf=(3/2*μ*v)^(1/2)
println("Feynman Polaron Radius: $rf , electron-mass unit")
println("Assuming Bohr, that's ",rf*5.2918e-11," m")
rf=(3/2*μ*v*MassElectron)^(1/2)
println("Feynman Polaron Radius: $rf, SI")
# Units?!? Bohr?
In [27]:
# ((freq THz)) ((IR Activity / e^2 amu^-1))
# These data from MAPbI3-Cubic_PeakTable.csv
# https://github.com/WMD-group/Phonons/tree/master/2015_MAPbI3/SimulatedSpectra
# Data published in Brivio2015 (PRB)
# https://doi.org/10.1103/PhysRevB.92.144308
MAPI= [
96.20813558773261 0.4996300522819191
93.13630357703363 1.7139631746083817
92.87834578121567 0.60108592692181
92.4847918585963 0.0058228799414729
92.26701437594754 0.100590086574602
89.43972834606603 0.006278895133832249
46.89209141511332 0.2460894564364346
46.420949316788 0.14174282581124137
44.0380222871706 0.1987196948553428
42.89702947649343 0.011159939465770681
42.67180170168193 0.02557751102757614
41.46971205834201 0.012555230726601503
37.08982543385215 0.00107488277468418
36.53555265689563 0.02126940080871224
30.20608114002676 0.009019481779712388
27.374810898415028 0.03994453721421388
26.363055017011728 0.05011922682554448
9.522966890022039 0.00075631870522737
4.016471586720514 0.08168931020200264
3.887605410774121 0.006311654262282101
3.5313112232401513 0.05353548710183397
2.755392921480459 0.021303020776321225
2.4380741812443247 0.23162784335484837
2.2490917637719408 0.2622203718355982
2.079632190634424 0.23382298607799906
2.0336707697261187 0.0623239656843172
1.5673011873879714 0.0367465760261409
1.0188379384951798 0.0126328938653956
1.0022960504442775 0.006817361620021601
0.9970130778462072 0.0103757951973341
0.9201781906386209 0.01095811116040592
0.800604081794174 0.0016830270365341532
0.5738689505255512 0.00646428491253749
#0.022939578929507105 8.355742795827834e-05 # Acoustic modes!
#0.04882611767873102 8.309858592685e-06
#0.07575149723846182 2.778248540373041e-05
]
# Change to SI, but not actually needed as units cancel everywhere
#MAPI_SI = [ MAPI_orig[:,1].*10^12*2*π MAPI_orig[:,2].*1 ]
# OK, black magic here - perhaps our units of oscillator strength are not what we need? maybe already effectively 'squared'?
#MAPI = [ MAPI[:,1] MAPI[:,2].^0.5]
MAPI_low=MAPI[19:33,:] # Just inorganic components, everything below 10THz; modes 3-18
Out[27]:
In [28]:
# Hellwarth Table II - BiSiO frequencies
HellwarthII = [
106.23 8.86
160.51 9.50
180.33 20.85
206.69 10.05
252.76 27.00
369.64 61.78
501.71 52.87
553.60 86.18
585.36 75.41
607.29 98.15
834.53 89.36
]
Out[28]:
In [29]:
# Most simple scheme
# Hellwarth (58), assuming further typo on LHS, actually should be W_e
function HellwarthBscheme(LO)
H58 = sum( (LO[:,2].^2)./ LO[:,1].^2 )
println("Hellwarth (58) summation: ",H58)
H59 = sum( LO[:,2].^2 ) # sum of total ir activity squarred
println("Hellwarth (59) summation (total ir activity ^2): ", H59)
println("Hellwarth (59) W_e (total ir activity ): ", sqrt(H59))
omega = sqrt(H59 / H58)
println("Hellwarth (61) Omega (freq): ",omega)
end
HellwarthBscheme(HellwarthII)
println(" should agree with values given in Hellwarth(60) W_e=196.9 cm^-1 and Hellwarth(61) Ω_e=500 cm^-1")
println("\t MAPI: (all values)")
HellwarthBscheme(MAPI)
println("\t MAPI: (low-frequency, non molecular IR)")
HellwarthBscheme(MAPI_low)
In [30]:
# Check this though, take the low-freq components and:
MAPI_low=[
9.522966890022039 0.00075631870522737
4.016471586720514 0.08168931020200264
3.887605410774121 0.006311654262282101
3.5313112232401513 0.05353548710183397
2.755392921480459 0.021303020776321225
2.4380741812443247 0.23162784335484837
2.2490917637719408 0.2622203718355982
2.079632190634424 0.23382298607799906
2.0336707697261187 0.0623239656843172
1.5673011873879714 0.0367465760261409
1.0188379384951798 0.0126328938653956
1.0022960504442775 0.006817361620021601
0.9970130778462072 0.0103757951973341
0.9201781906386209 0.01095811116040592
0.800604081794174 0.0016830270365341532
0.5738689505255512 0.00646428491253749
#0.022939578929507105 8.355742795827834e-05
#0.04882611767873102 8.309858592685e-06
#0.07575149723846182 2.778248540373041e-05
]
H58 = sum( (MAPI_low[:,2].^2)./ MAPI_low[:,1].^2 )
println("Hellwarth (58) summation: ",H58)
H59 = sum( MAPI_low[:,2].^2 ) # sum of total ir activity squarred
println("Hellwarth (59) summation (total ir activity ^2): ", H59)
omega = sqrt(H59 / H58)
println("Hellwarth (61) Omega (Thz): ",omega)
# That feels like it should be the right answer... for whatever reason the Hellwarth
# prescription isn't cancelling out the high-frequency components as you would expect
# it to from an integration across the Lorentz oscillators
In [ ]:
In [31]:
# More complex scheme, involving thermodynamic Beta
# Hellwarth(50), RHS
const hbar = const ħ = 1.05457162825e-34; # kg m2 / s
const eV = const q = const ElectronVolt = 1.602176487e-19; # kg m2 / s2
const me=MassElectron = 9.10938188e-31; # kg
const Boltzmann = const kB = 1.3806504e-23; # kg m2 / K s2
function HellwarthAscheme(LO,T=295)
β=LO[:,1].*2*pi*1E12*ħ/(kB*T) #assuming units of THz
H50 = sum( ((LO[:,2].^2).*coth(β))./LO[:,1] )
println("Hellwarth (50) summation: ",H50)
H51= sum( LO[:,2].^2 ) # sum of total ir activity squarred
println("Hellwarth (51) summation (total ir activity ^2): ", H51)
println("Hellwarth (51) W_e (total ir activity ): ", sqrt(H51))
# OK; so this is deriving Omega / coth(Beta/2)
omegacoth=H51/H50
println("omegacoth: ",omegacoth)
# NOT FINISHED - need to somehow decouple Omega from both sides of the eqn.
for freq in 0.25:0.25:20
pseudo_omega=omegacoth*coth(freq * 2*pi*1E12*ħ/(2*kB*T))
println("freq: $freq pseudo-omega: $pseudo_omega")
end
end
HellwarthAscheme(HellwarthII)
HellwarthAscheme( [HellwarthII[:,1].*0.02998 HellwarthII[:,2]] ) # convert data to Thz
println(" should agree with values given in Hellwarth\n TableII: H50sum= 91.34 cm^-1, \n W_e=196.9 cm^-1 and Hellwarth(53) Ω_e=504 cm^-1")
println("\t MAPI: (all values)")
HellwarthAscheme(MAPI)
println("\t MAPI: (low-frequency, non molecular IR)")
HellwarthAscheme(MAPI_low)
In [32]:
# Comparison via Hellwarth (37)
# MAPI 4.5, 24.1, 2.25THz - 75 cm^-1 ; α=
ε_Inf=4.5
ε_S=24.1
(1/ε_Inf - 1/ε_S) #/(4*pi)
# By Hellwarth (37) this should be equivalent to the summation Hellwarth(58) above
Out[32]:
In [33]:
# Hellwarth 44 defin. for alpha
Ry=13.605693*q
effectivemass=0.12
α=(1/ε_Inf - 1/ε_S) * sqrt(effectivemass * Ry / (ħ*2.25E12*2*pi))
Out[33]:
In [34]:
# Integrate through Lorentz oscillators to get dielectric fn
# Should give 'extra' contribution from these modes, extrapolated to zero omega
function integrate_dielectric(LO,V0)
summate=sum( (LO[:,2])./(LO[:,1].^2) )
summate*4*π/V0
end
Å=1E-10 # angstrom in metres
r=6.29Å # Sensible cubic cell size
V0=(r)^3
println("volume: $V0")
const amu=1.66054e-27
const ε0=8.854187817E-12
MAPI_SI = [ MAPI[:,1].*10^12*2*π MAPI[:,2]./(q^2/amu) ]
println(" MAPI: ",integrate_dielectric(MAPI,1.0))
println(" MAPI_low: ",integrate_dielectric(MAPI_low,1.0))
println(" MAPI_SI: ",integrate_dielectric(MAPI_SI,V0))
println(" MAPI_SI: fudged epislon0 ",integrate_dielectric(MAPI_SI,V0)*ε0/(4*π))
println(" MAPI_SI_low: fudged epislon0 ",integrate_dielectric(MAPI_SI[19:33,:],V0)*ε0/(4*π))
println()
println("From ε_S-ε_Inf, expect this to be: ",ε_S-ε_Inf)
In [35]:
sum(MAPI[:,2]./MAPI[:,1].^2 )
Out[35]:
In [36]:
sum(MAPI_low[:,2]./MAPI_low[:,1].^2 )
Out[36]:
In [41]:
println("Debye temperatures...")
freq=2.25E12 # Thz
ω = (2*pi)*freq
println("MAPI: 2.25 THz: ",hbar*ω/kB)
# GaP:
# 10.1002/pssb.2221200123
# TO mode at Gamma: 402.5 cm^-1
freq=12.06665E12 # Thz
ω = (2*pi)*freq
hbar*ω/kB
println("GaP: 12.066 THz: ",hbar*ω/kB)
In [56]:
meSmallAlpha(α)=α/6 + 0.025*α^2
# (46) In Feynman1955
meLargeAlpha(α)=16*α^4 / (81*π^4)
meLargeAlpha(α)=202*(α/10)^4
println("MAPI: electron, α=2.40 meSmallAlpha(α)=",meSmallAlpha(2.40)," meLargeAlpha(α)=",meLargeAlpha(2.40))
println("MAPI: hole, α=2.68 meSmallAlpha(α)=",meSmallAlpha(2.68)," meLargeAlpha(α)=",meLargeAlpha(2.68))
In [ ]:
In [ ]:
In [ ]: