TIDY3D
LEARNING CENTER

Lectures

Tutorials

Applications

Modeling dispersive material in FDTD

By Weiliang Jin, Zongfu Yu and Shanhui Fan

Material 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.

Share On:

Visit public projects:

Download .ipynb files: 1 2 View Presentation Slides

Additional information: This was updated in Apr 13, 2022

Tutorial 5: part 1 – surface plasmon polariton resonance

In this notebook, we show how to include simple analytical dispersion relations. As a detailed example, we will walk you through the simulation of surface plasmon polaritons (SPPs) at the gold-air interface, where the response of gold is approximated by the Drude model.

Setup

First we import Tidy3D and the other packages needed. If it is not installed, please perform the installation and account authentication by following the detailed instructions here.

# standard python imports
import matplotlib.pyplot as plt
import numpy as np
import tidy3d as td
from tidy3d import web
from tidy3d.constants import C_0, HBAR

td.config.local_cache.enabled = False

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=r"Re[$\epsilon$]")
ax.plot(lambda_list, ep.imag, label=r"Im[$\epsilon$]")
ax.set_xlabel(r"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

# dipole 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(r"x ($\mu$m)")
ax.set_ylabel(r"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="data/data_gold_spp.hdf5")

16:10:42 UTC Created task 'lecture05_gold_SPP' with resource_id                 
             'fdve-2c5f18fc-8525-48b5-9160-a69c8bbfbe0e' and task_type 'FDTD'.

             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-2c5f18fc-852
             5-48b5-9160-a69c8bbfbe0e'.

             Task folder: 'default'.

             Estimated FlexCredit cost: 0.063. This assumes the FDTD solver runs
             for the full simulation time; if early shutoff is reached, the     
             billed cost can be lower. Use 'web.real_cost(task_id)' to get the  
             billed FlexCredit cost after a simulation run.

             status = queued

             To cancel the simulation, use 'web.abort(task_id)' or              
             'web.delete(task_id)' or abort/delete the task in the web UI.      
             Terminating the Python script will not stop the job running on the 
             cloud.

16:10:52 UTC status = preprocess

16:10:56 UTC starting up solver

             running solver

16:11:11 UTC early shutoff detected at 69%, exiting.

             status = postprocess

16:11:13 UTC status = success

16:11:15 UTC View simulation result at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-2c5f18fc-852
             5-48b5-9160-a69c8bbfbe0e'.

             Loading results from 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(r"x ($\mu$m)")
ax[0].set_ylabel(r"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(r"x ($\mu$m)")
ax[1].set_ylabel(r"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.

Setup

First we import Tidy3D and the other packages needed. If it is not installed, please perform the installation and account authentication by following the detailed instructions here.

# standard python imports
import matplotlib.pyplot as plt
import numpy as np
import tidy3d as td
from tidy3d import web

# fitter plugins
from tidy3d.plugins.dispersion import FastDispersionFitter as Fitter

td.config.local_cache.enabled = False

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=database/data/main/Si/nk/Green-2008.yml"
)
fitter = fitter.copy(update={"wvl_range": (0.4, 1.2)})  # um

Let’s visualize 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(r"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(max_num_poles=2, tolerance_rms=5e-2)
medium = medium.copy(update={"name": "fitted silicon"})

Best weighted RMS error: 0.0392 ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 0:00:00

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

for count_slab, (t, mat) in enumerate(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

# 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(r"x ($\mu$m)")
ax.set_ylabel(r"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")

21:10:54 UTC Created task 'lecture05_silicon' with resource_id                  
             'fdve-843ca463-9e17-44a1-8145-ca3079792897' and task_type 'FDTD'.  
             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-843ca463-9e1
             7-44a1-8145-ca3079792897'.                                         
             Task folder: 'default'.                                            
↑ simulation.hdf5.gz ━━━━━━━━━━━━━━━━━━━━━━━━━ 100.0% • 4.6/4.6 kB • ? • 0:00:00

21:10:56 UTC Estimated FlexCredit cost: 0.154. This assumes the FDTD solver runs
             for the full simulation time; if early shutoff is reached, the     
             billed cost can be lower. Use 'web.real_cost(task_id)' to get the  
             billed FlexCredit cost after a simulation run.                     
21:11:00 UTC status = queued                                                    
             To cancel the simulation, use 'web.abort(task_id)' or              
             'web.delete(task_id)' or abort/delete the task in the web UI.      
             Terminating the Python script will not stop the job running on the 
             cloud.                                                             
21:12:02 UTC starting up solver                                                 
             running solver                                                     
21:12:29 UTC early shutoff detected at 72%, exiting.                            
solver progress (field decay = 1.30e-05) ━━━━━━━━━━━━━━━━━━━━━━━━━━ 100% 0:00:00

21:12:30 UTC status = success                                                   
             View simulation result at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-843ca463-9e1
             7-44a1-8145-ca3079792897'.                                         
↓ simulation_data.hdf5.gz ━━━━━━━━━━━━━━━━━━ 100.0% • 10.2/10.2 kB • ? • 0:00:00

             Loading results from 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(r"Wavelength ($\mu m$)")
ax.set_ylabel("Transmission")
ax.set_xlim(fitter.wvl_range)
ax.legend()
plt.show()

Download FDTD Modeling Material Dispersion

TIDY3D LEARNING CENTER

Modeling dispersive material in FDTD

Tutorial 5: part 1 – surface plasmon polariton resonance

Setup

Gold Dispersion with Drude Model

Other Basic Parameters

Create Simulation

Plot the Structure

Running Simulation

Post Run Analysis

Tutorial 5: part 2 – silicon slab transmission

Setup

Dispersion Fitter

Other Basic Setups

Create Simulation

Plot the Structure

Run the Simulation

Postprocess and Plot

Cookie preferences

Privacy choices

TIDY3D
LEARNING CENTER