DISPERSION FITTING TUTORIAL

NANOCPP allows to add materials with an arbitrary dispersion, thanks to a formulation based on an unconditionally stable algorithm with the same veocity of explicit schemes.

Below you will find a tutorial on how to fit your own dispersion model.

In [1]:
import os, sys
from matplotlib import pyplot as plt
import numpy as np
import math
from scipy.optimize import curve_fit
from scipy.signal import residue
import random as rnd
from importlib import reload
import warnings
warnings.filterwarnings('ignore')

The standard procedure involved in the simulation of dispersive materials requires a nonlinear analytical representation of the dielectric function ε (ω) expressed in terms of a set of known functions. Historically, three models were available: the Debye relaxation, the Lorentz resonance and the Drude model of metals. These models are defined in terms of a set of parameters related to the physical phenomena occurring in the medium (e.g., intra-band transitions, infrared absorption). In each model the dielectric function is defined as a sum of complex poles or poles-pairs, whose analytical representation mimics the relevant features of the dielectric function of the medium (e.g., relaxation dynamics, resonances).

Their expressions read

$$ \begin{equation} \varepsilon_{D e b y e}(\omega)=\varepsilon_{\infty}+\sum_{k=1}^{N} \frac{\Delta \varepsilon_{k}}{1+j \omega \tau_{k}} \end{equation} $$

$$ \begin{equation} \varepsilon_{\text {Lorent}}(\omega)=\varepsilon_{\infty}+\sum_{k=1}^{N} \frac{\Delta \varepsilon_{k} \omega_{p}^{2}}{\omega_{p}^{2}+2 j \omega \delta_{k}-\omega^{2}} \end{equation} $$

$$ \begin{equation} \varepsilon_{\text {Drude}}(\omega)=\varepsilon_{\infty}+\sum_{k=1}^{N} \frac{\omega_{k}^{2}}{\omega^{2}-j \omega \gamma_{k}} \end{equation} $$

These models, when considered singularly, are accurate only for a narrow frequency range for every set of parameters. In order to obtain a broadband model a sum of many terms is required in order to obtain an accurate representation of the dielectric function. More advanced techniques involve the introduction of mixed terms, as in the case of the Lorentz-Drude model, where Lorentz and Drude poles are coupled together, and whose applications are mainly in the analytical representation of metals. While these approaches provide a higher degree of accuracy, they usually require a large number of terms. In the case of an FDTD simulation, hence, a multi-pole formulation requires the computation and the storing of a large number of complex variables, significantly increasing the computational complexity

Recently, a new analytical method, known as the Critical Point (CP) model, has been proposed to overcome these difficulties.

$$ \begin{equation} \varepsilon_{C P}(\omega)=\varepsilon_{\infty}+\sum_{l=0}^{N}\left(\frac{c_{l}}{j \omega-d_{l}}+\frac{c_{l}^{*}}{j \omega-d_{l}^{*}}\right) \end{equation} $$

Iterating and simplifying, the dielectric function $\varepsilon_{CP}(\omega)$ becomes

$$ \begin{equation} \varepsilon_{C P}(\omega)=\frac{\sum b_{l}(j \omega)^{l}}{\sum a_{l}(j \omega)^{l}} \end{equation} $$

In NANOCPP, the dispersion model is expressed exactly in this form.

In [2]:
os.chdir('../dispersion/')
import materials as mt
In [3]:
def reload_fdtd():
    global tdcmt, myfun, frac_fit, exp_fit, disp_fit, mt, unpack, harm
    import fdtd_utils
    reload(fdtd_utils)
    import materials as mt
    import colors
    reload(mt)
    reload(colors)
    from fdtd_utils import tdcmt
    from fdtd_utils import myfun
    from fdtd_utils import frac_fit
    from fdtd_utils import exp_fit
    from fdtd_utils import disp_fit
    from fdtd_utils import unpack
    from fdtd_utils import harm

It is important to make sure that your model satisfy a few condition necessary for FDTD simulations to be stable.

For each simulation, it is important to calculate the critical frequency = 1/π/dt = 1/π/(dx·C/c), where dx is the resolution and c is the speed of light. The Courant factor C = 0.5 for a default 3D simulation. So for example, at resolution dx = 50 nm, we get the critical frequency in UV region: 1/3.14/(5e-8*.5/3e8) = 3.82e15 Hz.

Based on some experimenting with NANOCPP, it can be observed that following three conditions are needed for stability of the simulation:

  • At the critical frequency and above, the permittivity must be positive.

  • No Lorentzian oscillator may be defined above critical frequency, and preferably not above 1/2 of it to be safe.

  • (Make sure your geometry provides enough space between any resonator and perfectly matched layers.)

fdtd_utils allows you to construct a model using both local and remote data (from filmmetrics website).

In [4]:
reload_fdtd()
#si = mt.material('Si', source = 'web')
# **kargs takes collects the parameters using for loading data (for instance skiprows)
si = mt.material('Si.txt', source = 'file', skiprows=1)

All wavelengths needs to be in nanometers by default.

In [5]:
si.wl[::20]
Out[5]:
array([  210.01,   410.03,   610.04,   810.05,  1010.07])

For fitting you can use disp_fit class. The units parameters used below determines the output wavelengths as $\frac{\lambda}{units}$. In this particular example, the wavelength for the material it is measured in nm. By using units=1000 we assume that the output wavelength are in units of $\mu$m. If we want SI units we need to select $units=10^9$.

In [7]:
plt.figure(figsize=(5,2.5))
plt.xlabel('omega')
#plt.ylabel('$\Re$ $\varepsilon_r$, $\Im$ $\varepsilon_r$')
fit = disp_fit(si)
#see the notes below
fit.fit(N=1,units=1e3)

critical w (2pi/c) = 1.800633e-01, epsilon(3w/4) = (1.084178e+00,-9.043958e-02)

During the fit function call, the resulting string in NANOCPP format, corresponding to the polynomial coefficients of the calculated model, is automatically saved in the str.txt file.

In [8]:
!cat str.txt

3 5.834e+32 4.196e+22 2.072e+12 1.000e+00 5.242e+31 3.770e+21 1.884e+12 1.000e+00

Test

Poles

We here investigate the stability of the model based on the criteria described in the introduction. Namely, a sign of a real part of dielectric permittivity and poles position. We use the function check_FDTD and pass the grid resolution used in the FDTD simulations.

In [8]:
fit.check_FDTD(dx=5)

critical w (2pi/c) = 1.800633e-01, epsilon(3w/4) = (1.086907e+00,-5.150791e-02)

Let's plot $e_r, e_i$ vs. $\frac{w - w_{cr}}{w_{cr}}$ to see function behavior in the realm of critical omega

In [9]:
wcr = fit.wcr
w0 = fit.wn[0]
w1 = fit.wn[-1]
w = np.linspace(wcr,100*wcr,100)
plt.figure(figsize=(10,5), dpi = 200)
plt.xlabel(r'$\frac{w - w_{cr}}{w_{cr}}$')
fit.plot(w)
#value at critical omega
fit.out(wcr)
Out[9]:
(1.0883566426090829-0.045729021130875631j)

The position of Lorenzian pole can be calculated using the poles in residue expansion.

$$\begin{equation} \varepsilon(\omega)=\varepsilon_{\infty}+\sum_{l=0}^{N}\left(\frac{c_{l}}{j \omega-d_{l}}+\frac{c_{l}^{*}}{j \omega-d_{l}^{*}}\right) \end{equation}$$

$$\begin{equation} \varepsilon(\omega)=\varepsilon_{\infty}+\sum_{l=0}^{N}\left[\frac{2\left(\mathfrak{R}\left\{c_{l}\right\}-\mathfrak{R}\left\{d_{l}\right\}\right)-2 j \omega \Re\left\{c_{l}\right\}}{-\omega^{2}+\left|d_{l}\right|^{2}-2 j \omega \Re\left\{d_{l}\right\}}\right] \end{equation}$$

$$ \begin{equation} \varepsilon_{\text {Lornenz}}(\omega)=\varepsilon_{\infty}+\sum_{k=1}^{N} \frac{\Delta \varepsilon_{k} \omega_{p}^{2}}{\omega_{p}^{2}+2 j \omega \delta_{k}-\omega^{2}} \end{equation} $$

$$ \omega_p^2 = |d|^2$$

In [10]:
d, c, e_inf = residue(fit.a,fit.b)
np.array_split(d,2)
np.abs(d) / wcr
Out[10]:
array([ 0.0214469 ,  0.0214469 ,  0.07961505,  0.07961505,  0.52621233])

This completes the pre-analysis of the fitted model. Now we can start the simulation.

$ python test.py res_ex > output.txt

The test.py script will generate 3 folders:

  • res_ex
  • res_ex_cw
  • res_ex_harm

The first folder contains the results with an spl source for the dielectric slab structure.

The second folder contains the results with multiple spl sources for the dielectric slab.

Fitting

In [14]:
plt.figure(figsize=(10,5))
n = 3
resdir = 'res_ex/'
freqs = np.fromfile(resdir + 'epsilon_imag.bin')[::n+1]
wl = 2*np.pi*3e8 / freqs
si_re = np.fromfile(resdir + 'epsilon_real.bin')[2::n+1]
si_im = np.fromfile(resdir + 'epsilon_imag.bin')[2::n+1]
plt.xlim(0.2,1)
plt.xlabel('Wavelength, $\mu m$')
plt.plot(wl, si_re)
plt.plot(wl, si_im)
plt.plot(si.wl / 1e3, si.er, '-o')
plt.plot(si.wl / 1e3, -si.ei, '-o')
plt.show()

Energy distribution

It is useful to also check an energy distribution in the whole space, because sometimes instabilities could be localized.

In [18]:
src = np.fromfile(resdir_disp+'erg_src.bin')
src = src.reshape(161,161)
Exz = np.fromfile(resdir_disp+'erg_1.bin')
Exz = Exz.reshape(161,161)
plt.imshow(src)
Out[18]:
<matplotlib.image.AxesImage at 0x14ca8f3cd278>
In [19]:
Exz.max(), Exz.min()
Out[19]:
(1.5137370741555472e-24, 3.3716960903026535e-33)

The idea for the last test is to measure the actual material absorption (complex part of refractive index k) and (real part n) by checking a spatial distribution of electric field.

Real and complex components of refractive index expressed through corresponding components of dielectric permittivity as follows

$$ \begin{equation} n=\sqrt{\frac{\left|\varepsilon_{\mathrm{r}}\right|+\varepsilon^{'}}{2}} \end{equation} $$

$$ \begin{equation} \kappa=\sqrt{\frac{\left|\varepsilon_{\mathrm{r}}\right|-\varepsilon^{'}}{2}} \end{equation} $$

For this test we can consider a geometry made by a simple silicon slab.

In homogeneous material, a component of plane wave expressed as a sum of two complex conjugated exponents (or what is equivalent has ~cos(nkz) dependency).

$$ E = \frac{1}{2}\hat{E}\{ e^{inkz} + e^{-inkz} \} $$

In [23]:
ex = np.fromfile(resdir_disp_cw + 'ex-10.bin', dtype=np.float64)
ex = ex.reshape(161,402)
plt.imshow(ex)
plt.xlabel('z')
plt.ylabel('x')
plt.show()

Considering field distribution along z, and fitting it with a sum of exponential, it is possible to extract both real and imaginaru part of refractive index.

In [24]:
y = ex[80,230:-50]
plt.plot(y)
plt.show()

The fitting procedure will be done using exp_fit function from fdtd_utils.py. Below you can find an implementation of this procedure.

In [25]:
wl2 = np.linspace(0.2,1.0,20)
def get_n_k(ex):
    y = ex[80,230:-50]
    x=np.linspace(0,1.,len(y))
    y1=y[:]
    x1=x[:]
    b = exp_fit()
    b.mtm(y1,2,x1[2]-x1[1])
    return b.alfa[0], b.w[0]
In [26]:
res = []
for i in range(20):
    ex = np.fromfile(resdir_disp_cw + 'ex-%i.bin' % i, dtype=np.float64)
    ex = ex.reshape(161,402)
    n, k = get_n_k(ex)
    res.append([n, k])
res = np.array(res)
In [27]:
n = 3
freqs = np.fromfile(resdir + 'epsilon_imag.bin')[::n+1]
freqs = np.fromfile(resdir + 'epsilon_imag.bin')[::n+1]
wl = 2*np.pi*3e8 / freqs
si_re = np.fromfile(resdir + 'epsilon_real.bin')[2::n+1]
si_im = np.fromfile(resdir + 'epsilon_imag.bin')[2::n+1]
wl = 2*np.pi*3e8 / freqs
e_r = si_re + 1j*si_im
n = np.sqrt((np.abs(e_r)+np.real(e_r))/2.)
k = np.sqrt((np.abs(e_r)-np.real(e_r))/2.)
In [28]:
plt.xlim(0.2,1)
const = 3.8
plt.plot(wl2,-res[:,0] * wl2 / const,'ro')
plt.plot(wl2,res[:,1] * wl2 / const,'bo')
plt.plot(wl,n)
plt.plot(wl,k)
plt.xlabel('wavelength, um')
plt.show()