The complex dielectric function $\epsilon$ is given by $\epsilon = \epsilon_\infty + \sum \frac{f_0\omega_0^2}{\omega_0^2 - \omega^2 + i\gamma \omega}$
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate
from scipy.interpolate import interp1d
%matplotlib ipympl
################# simulation constants #######################
Einf = 1
w0 = [1000,1200]
g = [100,1020]
wStart = 100
wEnd = 10000
wStep = 1
f0 = [-1,-10]
wRange = np.arange(wStart,wEnd,wStep)
E = np.zeros(len(wRange))
def getLorenzian(w0,g,f0,wRange):
temp =f0*w0**2/(w0**2 - np.square(wRange) + 1j*g*wRange)
return [np.real(temp),np.imag(temp)]
for i in range(len(w0)):
E = E + getLorenzian(w0[i],g[i],f0[i],wRange)
E = E + Einf
plt.plot(wRange,E[0],label = 'Real part')
plt.plot(wRange,E[1],label = 'Imag part')
plt.legend()
plt.show()
Now, there has been numerous articles which talk about extracting n and k directly from reflectivity data. See (Determination of Optical Constants from Reflectance or Transmittance Measurements on Bulk Crystals or Thin Films,Evaluation of the One-Angle Reflection Technique for the Determination of Optical Constants)
The basic idea is as follows: Reflection is a complex quantity $\textbf{r} = re^{i\phi}$ . When the reflectivity is measured, we measure the absolute value of a complex quantity $R = |re^{i\phi}|^2$. R is measured as a function of wavelength/frequency of the incident radiation either in a UV-Vis spectrometer or similar equipment. Therefore, the complex reflection is just $\textbf{r} = \sqrt{R}e^{i\phi}$.
The complex refractive index $\textbf{n} = n - ik$ is related to complex reflectivity in this fashion \begin{equation} \textbf{r} = \frac{\textbf{n} - 1}{\textbf{n} + 1} \end{equation} in the case of normal incidence
So, individual components of complex relectivity is given by: \begin{align} n = \frac{1 - r^2}{ 1 - 2r cos\phi + r^2}\\ k = \frac{-2r sin\phi}{ 1 - 2r cos\phi + r^2}\\ \end{align}
So, if we measure experimentally R, then $r(w) = \sqrt{R(w)}$ and the phase \begin{equation} \phi(\omega) = \frac{-2\omega}{\pi} P\int_0^\infty \frac{ln\sqrt{R(w)}}{x^2 - \omega^2}dx \end{equation}
sData = np.loadtxt('SiReflectivity.csv',skiprows = 1,delimiter = ',')
fig1,ax1 = plt.subplots()
ax1.plot(sData[:,0],sData[:,5])
ax1.set_ylabel('Reflectivity')
Text(0, 0.5, 'Reflectivity')
rData = np.loadtxt('reflectivity MoSe2_Ga D1.txt',skiprows = 1)
w = 3E8/(rData[:,0]*1E-9) # 1/s
r = rData[:,1]/100
Now definiting a function which will take the reflectivity data and get me the n and k values
Defining calculation functions
def getTheta(rin,win):
phi = np.zeros(len(win))
phi_err = np.zeros(len(win))
rFun = interp1d(win,rin,bounds_error = False, fill_value=(rin[0],rin[-1]))
#Integrand = lambda w,wp: np.log(rFun(w))/(w + wp)
Integrand = lambda w,wp: np.log(np.sqrt(rFun(w)))/(w + wp)
for i in range(len(win)):
phi[i],phi_err = integrate.quad(Integrand,0,10, args = win[i],weight = 'cauchy',wvar = win[i],limit = 100,epsabs = 1E-20)#*2*win[i]/np.pi
#print(phi_err)
return phi*2*win/np.pi
def getnk(rin,theta):
n = np.zeros(len(rin))
k = np.zeros(len(rin))
for i in range(len(rin)):
#n[i] = (1 - rin[i]**2)/(1 - 2*rin[i]*np.cos(phi[i]) + rin[i]**2)
#k[i] = (-2*rin[i]*np.sin(phi[i]))/(1 - 2*rin[i]*np.cos(phi[i]) + rin[i]**2)
n[i] = (1 - rin[i])/(1 - 2*np.sqrt(rin[i])*np.cos(phi[i]) + rin[i])
k[i] = (-2*np.sqrt(rin[i])*np.sin(phi[i]))/(1 - 2*np.sqrt(rin[i])*np.cos(phi[i]) + rin[i])
return [n,k]
We will first calculate the n and k for Si to check our codes:
w = 3/(sData[:,0]) ## 10^14 1/s
r = sData[:,5]
phi = getTheta(r,w)
[nCalc,kCalc] = getnk(r,phi)
/tmp/ipykernel_761962/1633332591.py:8: IntegrationWarning: The maximum number of subdivisions (100) has been achieved. If increasing the limit yields no improvement it is advised to analyze the integrand in order to determine the difficulties. If the position of a local difficulty can be determined (singularity, discontinuity) one will probably gain from splitting up the interval and calling the integrator on the subranges. Perhaps a special-purpose integrator should be used. phi[i],phi_err = integrate.quad(Integrand,0,10, args = win[i],weight = 'cauchy',wvar = win[i],limit = 100,epsabs = 1E-20)#*2*win[i]/np.pi
## Lets plot the complex part of the reflectivity $\phi$
fig,ax = plt.subplots()
ax.plot(w,phi)
[<matplotlib.lines.Line2D at 0x7fc3106af2b0>]
Plotting the calculated values
fig2,ax2 = plt.subplots()
ax2.plot(sData[:,0],kCalc,label = 'k-Calculated')
ax2.plot(sData[:,0],nCalc,label = 'n-Calculated')
ax2.plot(sData[:,0],sData[:,4],label = 'k-reported')
ax2.plot(sData[:,0],sData[:,3],label = 'n-reported')
ax2.legend()
#ax1.set_ylim([-500,500])
<matplotlib.legend.Legend at 0x7fc310691d60>
inData = np.loadtxt('reflectivity MoSe2_Ga D1.txt',skiprows = 1)
w = 3/(inData[:,0]) # 1/s
r = inData[:,1]/100
phi = getTheta(r,w)
[nCalc,kCalc] = getnk(r,phi)
/tmp/ipykernel_761962/1633332591.py:8: IntegrationWarning: The maximum number of subdivisions (100) has been achieved. If increasing the limit yields no improvement it is advised to analyze the integrand in order to determine the difficulties. If the position of a local difficulty can be determined (singularity, discontinuity) one will probably gain from splitting up the interval and calling the integrator on the subranges. Perhaps a special-purpose integrator should be used. phi[i],phi_err = integrate.quad(Integrand,0,10, args = win[i],weight = 'cauchy',wvar = win[i],limit = 100,epsabs = 1E-20)#*2*win[i]/np.pi
fig3,ax3 = plt.subplots()
ax3.plot(inData[:,0],-kCalc,label = 'k')
ax3.plot(inData[:,0],nCalc,label = 'n')
ax3.set_xlabel('wavelength \(nm\)')
ax3.set_ylabel('complex refractive index')
ax3.legend()
<matplotlib.legend.Legend at 0x7fc3104f7dc0>
fig4,ax4 = plt.subplots()
ax4.plot(inData[:,0],-4*np.pi*kCalc*1E9/inData[:,0],label = 'Absorption Coefficient')
ax4.set_title('Absorption Coefficient')
ax4.set_xlim([400,900])
(400.0, 900.0)