Modeling dispersive material in FDTD
By Weiliang Jin, Zongfu Yu and Shanhui FanMaterial dispersion is a very common phenomenon in which the material responds differently to light of different color. In this lecture, we show how to include material dispersion in FDTD simulations.
- Show how to include simple analytical dispersion relations in FDTD, illustrated by an example of simulating surface plasmon polaritonic (SPP) resonance between gold and air interface.
- Introduction to a popular method for describing material dispersion in FDTD, known as the complex-conjugate pole-residue method.
- Show how to include more complicated dispersion relations that need to be inferred from tabulated data, illustrated by an example of simulating crystalline silicon slab transmission with the help of dispersion fitting tools.
Tutorial 5: Part 1 -- surface plasmon polariton resonance¶
In this notebook, we show how to include simple analytical dispersion relations. As an detailed example, we will walk you through the simulation of surface plasmon polaritons (SPPs) at the gold-air interface, where the reponse of gold is approximated by the Drude model.
# standard python imports
import numpy as np
import matplotlib.pyplot as plt
import tidy3d as td
from tidy3d import web
from tidy3d.constants import C_0, HBAR
Gold dispersion with Drude model¶
In the visible frequency range, the dispersion of gold can be well approximated by the Drude model.
# Drude parameters for gold
eps_inf = 9.84 # relative permittivity at infinite frequency
wp = 9.01 #eV, plasma frequency
gamma = 0.072 #eV, damping rate
# Unit conversion from ev to Hz
wp_hz = wp / HBAR / 2 / np.pi
gamma_hz = gamma / HBAR / 2 / np.pi
# Tidy3d Drude model
mat_gold = td.Drude(eps_inf = eps_inf, coeffs = [(wp_hz,gamma_hz)], name='gold Drude')
Visualize the dispersion
# Compute permittivity in a range of wavelengths
lambda_list = np.linspace(0.3,0.6,500) # um
freq_list = C_0/lambda_list
ep = mat_gold.eps_model(freq_list)
# Visualize the real and imaginary part of permittivity
plt.rcParams.update({'font.size': 16})
fig, ax = plt.subplots(1,tight_layout=True,figsize=(5,4))
ax.plot(lambda_list,ep.real,label="Re[$\epsilon$]")
ax.plot(lambda_list,ep.imag,label="Im[$\epsilon$]")
ax.set_xlabel('Wavelength ($\mu$m)')
ax.set_xlim(lambda_list[0],lambda_list[-1])
ax.legend()
plt.show()
Other basic parameters¶
Now let us define the rest of basic parameters to setup the simulation. We consider placing a dipole source near the gold-air interface to excite the SPP.
# thickness of gold slab and air space
t_gold = 0.5 #um
t_air = 4 #um
# diple source position
distance = 0.02 #um, distance between the dipole source and the gold-air interface
px = 1 # um, distance between the dipole source and the left boundary of the simulation domain
# simulation size
Lx = 10 #um. we'll be looking at SPP propagating along x-direction, so let's make it long along x-axis
Ly = 2 #um
Lz = t_gold+t_air
sim_size = Lx, Ly, Lz
# Wavelength and frequency range of the source
lambda_range = (0.4,0.6)
lambda0 = (lambda_range[0] + lambda_range[1])/2
freq0 = C_0/lambda0
freqw = 0.3*(C_0/lambda_range[0]-C_0/lambda_range[1])
# runtime
runtime = 10
t_stop = runtime/freqw
# frequencies and wavelengths of field monitor
# we consider two representative wavelengths: 0.45 um for mirror-like response, and 0.55 um for SPP response
monitor_lambdas = np.array([0.45,0.55]) #um
monitor_freqs = C_0 / monitor_lambdas
Create Simulation¶
Now we set everything else up (structures, sources, monitors, simulation) to run the example.
# gold slab
gold_slab = td.Structure(
geometry=td.Box(
center=(0, 0, -Lz),
size=(td.inf,td.inf,Lz+2*t_gold),
),
medium=mat_gold,
name='gold slab',
)
# dipole source
source = td.PointDipole(
source_time = td.GaussianPulse(
freq0=freq0,
fwidth=freqw
),
center=(-Lx/2+px, 0, -Lz/2+t_gold+distance),
polarization='Ex',
name='dipole',
)
# field monitor
monitor = td.FieldMonitor(
center = (0, 0, 0),
size = (Lx, 0, Lz),
freqs = np.sort(monitor_freqs),
name='field',
)
# nonuniform mesh
grid_spec=td.GridSpec.auto(
min_steps_per_wvl=20,
)
sim = td.Simulation(
center = (0, 0, 0),
size = sim_size,
grid_spec = grid_spec,
structures = [gold_slab],
sources = [source],
monitors = [monitor],
run_time = t_stop,
boundary_spec = td.BoundarySpec.all_sides(boundary=td.PML()),
)
Plot The Structure¶
Let's now plot the permittivity profile to confirm that the structure was defined correctly. We use the Simulation.plot() method to plot the materials only, which assigns a different color to each slab without knowledge of the material properties.
fig, ax = plt.subplots(1, tight_layout=True, figsize=(6, 4))
sim.plot(y=0, freq=freq0, ax=ax);
ax.set_title('')
ax.set_xlabel('x ($\mu$m)')
ax.set_ylabel('z ($\mu$m)')
plt.show()
Running simulation¶
We will submit the simulation to run as a new project.
sim_data = web.run(sim, task_name='lecture05_gold_SPP', path=f'data/data_gold_spp.hdf5')
Post Run Analysis¶
We visualize the field at the wavelength 450nm and 550 nm.
fig, ax = plt.subplots(1,2,figsize=(16, 4), tight_layout=True)
sim_data.plot_field('field', 'Ex', y=0, f=C_0/0.45, val='real', ax = ax[0], vmax = None)
ax[0].set_title('Wavelength: 450 nm')
ax[0].set_xlabel('x ($\mu$m)')
ax[0].set_ylabel('z ($\mu$m)')
sim_data.plot_field('field', 'Ex', y=0, f=C_0/0.55, val='real', ax = ax[1], vmax = None)
ax[1].set_title('Wavelength: 550 nm')
ax[1].set_xlabel('x ($\mu$m)')
ax[1].set_ylabel('z ($\mu$m)')
plt.show()
Tutorial 5: Part 2 -- silicon slab transmission¶
In this part, we move on to treat more complex material dispersion. As an example, we illustrate how to use our dispersion fitting tool to fit the dispersion relation of crystalline silicon based on tabulated data, and subsequently convert it to Tidy3D components for use in computing silicon slab transmission.
# standard python imports
import numpy as np
import matplotlib.pyplot as plt
import tidy3d as td
from tidy3d import web
# fitter plugins
from tidy3d.plugins.dispersion import StableDispersionFitter as Fitter, AdvancedFitterParam
Dispersion fitter¶
We load the dispersion data for crystalline silicon at room temperature directly from a URL. Considering that we are interested in the solar cell applications, we narrow the wavelength range to [0.4 um,1.2 um] to perform the fitting.
fitter = Fitter.from_url('https://refractiveindex.info/data_csv.php?datafile=data/main/Si/Green-2008.yml')
fitter = fitter.copy(update={'wvl_range' : (0.4, 1.2)}) #um
Let's visulize the complete dispersion data.
plt.rcParams.update({'font.size': 18})
fig, ax = plt.subplots(1,figsize=(5,4))
ax.scatter(fitter.wvl_um, fitter.n_data, label='n', color='crimson', edgecolors='black', linewidth=0.5)
ax.scatter(fitter.wvl_um, fitter.k_data, label='k', color='blueviolet', edgecolors='black', linewidth=0.5)
ax.set_xlabel('Wavelength ($\mu m$)')
ax.legend()
plt.show()
The fitting is then performed. More details on setting up the fitter can be found here.
medium, rms_error = fitter.fit(
num_poles=2,
tolerance_rms=5e-2,
num_tries=50,
)
medium = medium.copy(update={'name':"fitted silicon"})
Let’s visualize how well the fit captures the dispersion data in the wavelength range of interest. As you can see, this fit looks great and should be sufficient for our simulation.
fig, ax = plt.subplots(1,figsize=(6,5))
fitter.plot(medium,ax=ax)
plt.show()
Other basic setups¶
Now let us setup the rest of the basic components.
# Wavelength and frequency range
lambda_range = (0.4, 1.2)
lam0 = np.sum(lambda_range)/2
freq_range = (td.constants.C_0/lambda_range[1], td.constants.C_0/lambda_range[0])
# frequencies and wavelengths of monitor
Nfreq = 333
monitor_freqs = np.linspace(freq_range[0], freq_range[1], Nfreq)
monitor_lambdas = td.constants.C_0 / monitor_freqs
# central frequency, frequency pulse width and total running time
freq0 = monitor_freqs[Nfreq // 2]
freqw = 0.3 * (freq_range[1] - freq_range[0])
t_stop = 100 / freq0
# Thicknesses of slabs
t_slabs = [1.0] # um
# medium of slabs
mat_slabs = [medium]
# space between slabs and sources and PML
spacing = 1 * lambda_range[-1]
# simulation size
sim_size = Lx, Ly, Lz = (1, 1, 4*spacing + sum(t_slabs))
Create Simulation¶
Now we set everything else up (structures, sources, monitors, simulation) to run the example.
# First, we define the multilayer stack structure (a single layer in this example).
slabs = []
slab_position = -Lz/2 + 2*spacing
count_slab = 0
for t, mat in zip(t_slabs, mat_slabs):
slab = td.Structure(
geometry=td.Box(
center=(0, 0, slab_position + t/2),
size=(td.inf, td.inf, t),
),
medium=mat,
name='slab'+ str(count_slab),
)
slabs.append(slab)
slab_position += t
count_slab += 1
# Here we define the planewave source, placed just in advance (towards negative z) of the slab
source = td.PlaneWave(
source_time = td.GaussianPulse(
freq0=freq0,
fwidth=freqw
),
size=(td.inf, td.inf, 0),
center=(0, 0, -Lz/2+spacing*1.5),
direction='+',
pol_angle=0,
name='planewave',
)
# We are interested in measuring the transmitted flux, so we set it to be an oversized plane.
monitor = td.FluxMonitor(
center = (0, 0, Lz/2 - spacing),
size = (td.inf, td.inf, 0),
freqs = monitor_freqs,
name='trans',
)
Now it is time to define the simulation object.
sim = td.Simulation(
center = (0, 0, 0),
size = sim_size,
grid_spec=td.GridSpec.auto(min_steps_per_wvl=30),
structures = slabs,
sources = [source],
monitors = [monitor],
run_time = t_stop,
boundary_spec = td.BoundarySpec.pml(z=True),
)
Plot The Structure¶
plt.rcParams.update({'font.size': 12})
fig, ax = plt.subplots(1, tight_layout=True, figsize=(4, 6))
sim.plot(y=0, freq=freq0, ax=ax);
ax.set_title('')
ax.set_xlabel('x ($\mu$m)')
ax.set_ylabel('z ($\mu$m)')
plt.show()
Run the simulation¶
We will submit the simulation to run as a new project.
sim_data = web.run(sim, task_name='lecture05_silicon', path='data/data_silicon.hdf5')
Postprocess and Plot¶
Once the simulation has completed, we can download the results and load them into the simulation object.
Now, we compare the transmitted flux computed from a transfer matrix method (TMM) code to that computed from the FDTD simulation.
# import TMM package
import tmm
# prepare list of thicknesses including air boundaries
d_list = [np.inf] + t_slabs + [np.inf]
# loop through wavelength and record TMM computed transmission
transmission_tmm = []
for i, lam in enumerate(monitor_lambdas):
# create list of refractive index at this wavelength including outer material (air)
n_list = [1, np.sqrt(medium.eps_model(td.C_0/lam)), 1]
# get transmission at normal incidence
T = tmm.coh_tmm('s', n_list, d_list, 0, lam)['T']
transmission_tmm.append(T)
transmission = sim_data['trans'].flux
fig, ax = plt.subplots(1,figsize=(5,4))
ax.plot(monitor_lambdas, transmission_tmm, '-', label='TMM')
ax.plot(monitor_lambdas, transmission, 'k--', label='Tidy3D')
ax.set_xlabel('Wavelength ($\mu m$)')
ax.set_ylabel('Transmission')
ax.set_xlim(fitter.wvl_range)
ax.legend()
plt.show()
FDTD 101: Lecture 5
Today we will continue our tutorial. I'm Shanhui Fan from Flexcompute.
How to model dispersive materials in FDTD ?
We would talk about the modeling of dispersive material in FDTD method. Material dispersion is quite common, and this shows up when the dielectric permittivity or the dielectric function is a function of frequency. This shows up when the material responds differently to light at different colors or different frequencies. One common example of a dispersive material is metal in the optical frequency range. For example, the permittivity of gold can be very well described by the Drude model. The real part of the permittivity varies from being positive to negative in the wavelength range across the visible frequency spectrum. Dispersion of materials gives rise to interesting new physics that does not occur in non-dispersive materials.
Dispersion of gold gives rise to plasmonic polariton.
For example, in the case of gold this Drude model response gives rise to a strong plasmonic response. As a simple illustration of it, we can do a numerical experiment where we put a dipole source in the immediate vicinity of the interface between gold and vacuum. And we look at the field distribution on the slice as indicated by the gray plant here. When we choose a wavelength of 450 nanometer, you can see that the field radiates into vacuum because the gold essentially behaves like a mirror and that's perhaps what you might expect from your everyday experience of metal. However, if you lower the frequency or increase the wavelength a little bit to 550 nanometer. Now in this case, in addition to the radiation into vacuum. The dipole source here also excites a surface mode at the gold air interface. This surface mode is commonly referred to as the surface plasmon polariton and there's a very large literature about its physics and device application. I'm not going to review the physics associated with plasmonics. In fact, the study of these kinds of waves has given rise to the field of plasmonics.
What is the relationship between E field and D field in dispersive materials?
Up to this point in all the previous tutorials when we model materials, we assume the material to be dispersionless. For these materials, when you go from the electric field E at a particular time to the displacement field D, which contains the dipole moment of the material at the same time. What you do is that you simply multiply the electric field by the dielectric primitivity. In the time domain, the displacement field reacts instantaneously to the applied electric field. On the other hand, the dispersion of physics is somewhat different. I've shown you a descriptive description of the permittivity as a function of frequency. And this is a frequency domain description. Since we are doing the simulation in the time domain in the FDTD method, it would be useful to consider the corresponding time domain physics. There the multiplication of two functions in the frequency domain translates into a convolution operation in time domain. The displacement field is related to the electric field as a function of time convolved with a kernel. Dispersive material has a memory. Its displacement field depends not just on the electric field at current time, but also depends on the electric fields in the past. In the time domain modeling of a dispersive material, one needs to be able to account for this memory effect.
How to use complex conjugate pole residue methods to simulate dispersion in FDTD?
One of the popular methods is the complex conjugate pole residue pair method. In this method, you take the frequency permittivity and express them in this functional form consisting of a sum of poles. In these cases, for example, a_m is the pole that describes some sort of resonant frequencies of the material. As an example, for the Drude model that we have previously used to describe gold, you can use two pairs in this conjugate pole method to describe the Drude model. Next we will go a step deeper and look at the algorithm once you have described the permittivity in terms of these pole residuel pairs.
How to construct FDTD algorithms to describe the dispersion model?
For each of these poles, you define an auxiliary current to satisfy a differential equation that corresponds to this particular pole. This differential equation relates the auxiliary current to the applied electric field. Once you define all these auxiliary currents, you can then sum all their contributions together to describe the total dielectric response of the material. As you can see from this discussion, the computational cost here increases with the number of poles. In the example for a Drude model for metal, the Drude model naturally decomposes into these kinds of pole residue pairs and then one can use the algorithm as outlined here to describe it.
Simulate the transmission through a dispersive silicon slab in FDTD
In many practical situations you were typically given a tabulated data of the permittivity as a function of frequency or wavelengths. For example, you may be interested in modeling the behavior of a silicon slab in the visible and near infrared wavelength range. And in this case, you can go to standard references and they would tell you how the refractive index varies as a functional wavelength and these are not given as a functional form, but rather they are given in tabulated data. In this case, you need to use the pole residue pairs method to fit these tabulated data.
How to construct the dispersion model for silicon material in FDTD?
For silicon, if you are interested in modeling the behavior in the wavelength range spanning from 400 nanometer to 1.2 micron, which is a wavelength range approximately above the silicon band gap and so where the silicon is absorbing. Now in this case, as it turned out, only two poles are sufficient to give a very good modeling. So in Tidy3D you can find an utility that allows you to do fitting and here the dots here are the data taken from standard reference determined experimentally. And the solid line here is the fit. As we mentioned, this fit is important. And in fact, it's necessary in order to use the FDTD method to describe dispersive materials. Using a relatively small number of poles, it is possible to provide a fairly accurate description of most material dispersion over a fairly wide range of frequencies.
Comparing FDTD results with those from transfer matrix method for a dispersive silicon slab.
Once you have the capability to describe the dispersive material, it is possible to perform this broadband computation for dispersive structure in a single simulation. Here is an illustration where we compute the transmission spectrum of light normally incident upon a silicon slap over the wavelength range from 400 nanometer to 1.2 micron. We compare the result with a transfer matrix modeling, which is a frequency domain method. You can see the two agree very well. So to briefly summarize, I hope to give you an impression of how dispersive material is being modeled in the FDTD simulation.
