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

In [1]:

# 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 FastDispersionFitter as Fitter

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.

In [2]:

fitter = Fitter.from_url('https://refractiveindex.info/data_csv.php?datafile=database/data-nk/main/Si/Green-2008.yml')
fitter = fitter.copy(update={'wvl_range' : (0.4, 1.2)}) #um

Let's visulize the complete dispersion data.

In [3]:

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()

No description has been provided for this image

The fitting is then performed. More details on setting up the fitter can be found here.

In [4]:

medium, rms_error = fitter.fit(max_num_poles=2, tolerance_rms=5e-2)
medium = medium.copy(update={'name':"fitted silicon"})

Output()

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.

In [5]:

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.

In [6]:

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

In [7]:

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

In [8]:

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¶

In [9]:

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.

In [10]:

sim_data = web.run(sim, task_name='lecture05_silicon', path='data/data_silicon.hdf5')

17:24:06 PST Created task 'lecture05_silicon' with task_id
             'fdve-2e9adbd0-f29f-4e73-869d-e07360b8963e' and task_type 'FDTD'.

             View task using web UI at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-2e9adbd0-f29
             f-4e73-869d-e07360b8963e'.

Output()

17:24:07 PST status = queued

Output()

17:24:10 PST status = preprocess

17:24:14 PST Maximum FlexCredit cost: 0.156. Use 'web.real_cost(task_id)' to get
             the billed FlexCredit cost after a simulation run.

             starting up solver

17:24:15 PST running solver

             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.

Output()

17:24:38 PST early shutoff detected at 76%, exiting.

             status = postprocess

Output()

17:24:41 PST status = success

             View simulation result at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-2e9adbd0-f29
             f-4e73-869d-e07360b8963e'.

Output()

17:24:42 PST loading simulation 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.

In [11]:

# 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)

In [12]:

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()

In [ ]:

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¶

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¶

TIDY3D
LEARNING CENTER