Dispersive materials

Introduction / Setup

Here we show to to model dispersive materials in Tidy3D with an example showing transmission spectrum of a multilayer stack of slabs.

See on github, run on colab, or just follow along with the output below.

# get the most recent version of tidy3d
!pip install -q --upgrade tidy3d

# make sure notebook plots inline
%matplotlib inline

# standard python imports
import numpy as np
import matplotlib.pyplot as plt

import tidy3d as td
from tidy3d import web
Using Tidy3D credentials from stored file

First, let us define some basic parameters.

# Wavelength and frequency range
lambda_range = (0.5, 1.5)
lam0 = np.sum(lambda_range)/2
freq_range = (td.constants.C_0/lambda_range[1], td.constants.C_0/lambda_range[0])
Nfreq = 333

# frequencies and wavelengths of monitor
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 = [0.5, 0.2, 0.4, 0.3] # um

# Grid resolution (cells per um)
res = 150

# space between slabs and sources and PML
spacing = 1 * lambda_range[-1]

# simulation size
sim_size = Lx, Ly, Lz = (1.0, 1.0, 4*spacing + sum(t_slabs))

Defining Materials (4 Ways)

We now create materials for each slab in four differt ways. The first three of those show various ways that the Tidy3D Medium class can be used to define custom materials. The last one is a direct import from our material library.

  1. Simple, lossless dielectric defined by a real, dispersionless refracive index (or permittivity).

  2. Lossy material defined by real and imaginary part of the refractive index (\(n\)) and (\(k\)) at a given frequency or wavelength. Values are exact only at that frequency, so only good for narrow-band simulations.

  3. Simple, lossless dispersive material (one-pole fitting) defined by the real part of the refractive index \(n\) and the dispersion \(\mathrm{d}n/\mathrm{d}\lambda\) at a given frequency or wavelength. The dispersion must be negative.

  4. Dispersive material imported from our pre-defined library of materials.

More complicated dispersive materials can also be defined through their dispersive models (supported: Lorentz, Sellmeier, or Debye models). See here for more details on defining materials this way if the model parameters are known.

# simple, lossless, dispersionless material (either epsilon or n)
mat1 = td.Medium(epsilon=4.0)

# lossy material with n & k values at a specified frequency or wavelength
mat2 = td.Medium(n=3.0, k=0.1, freq=freq0)

# lossless dispersive material with n & dn/dlambda at a specified wavelength
mat3 = td.Sellmeier.from_dn(n=2.0, wl=lam0, dn_dwl=-0.02)

# dispersive material from tidy3d library
mat4 = td.material_library.BK7()

# put all together
mat_slabs = [mat1, mat2, mat3, mat4]

Create Simulation

Now we set everything else up (structures, sources, monitors, simulation) to run the example.

First, we define the multilayer stack structure.

slabs = []
slab_position = -Lz/2 + 2*spacing
for t, mat in zip(t_slabs, mat_slabs):
    slab = td.Box(
        center=(0, 0, slab_position + t/2),
        size=(td.inf, td.inf, t),
    slab_position += t

We must now define the excitation conditions and field monitors. We will excite the slab using a normally incident (along z) planewave, polarized along the x direciton.

# Here we define the planewave source, placed just in advance (towards negative z) of the slab
source = td.PlaneWave(
    source_time = td.GaussianPulse(
        frequency = freq0,
        fwidth = freqw
    position=-Lz/2 + spacing,

Here we define the field monitor, placed just past (towards positive z) of the stack.

# We are interested in measuring the transmitted flux, so we set it to be an oversized plane.
monitor = td.FreqMonitor(
    center = (0, 0, Lz/2 - spacing),
    size = (td.inf, td.inf, 0),
    freqs = monitor_freqs,

Now it is time to define the simulation object.

sim = td.Simulation(
    center = (0, 0, 0),
    size = sim_size,
    mesh_step = 1/res,
    structures = slabs,
    sources = [source],
    monitors = [monitor],
    run_time = t_stop,
    pml_layers = (0, 0, 15)
Initializing simulation...
Mesh step (micron): [6.67e-03, 6.67e-03, 6.67e-03].
Simulation domain in number of grid points: [150, 150, 1140].
Total number of computational grid points: 2.56e+07.
Total number of time steps: 21651.
Estimated data size (GB) of monitor monitor: 0.0000.
WARNING: Simulation frequency range exceeds the range of validity of the material model of material BK7.

Plot The Structure

Let’s now plot the permittivity profile to confirm that the structure was defined correctly.

First we use viz_mat_2D to plot the materials only, which assigns a different color to each slab without knowledge of the material properties.


Next, we use viz_eps_2D to vizualize the permittivity of the stack. However, because the stack contains dispersive materials, we need to specify the frequency of interest as an argument to the plotting tool. Here we show the permittivity at the lowest and highest frequences in the range of interest. Note that in this case, the real part of the permittivity (being plotted) only changes slightly between the two frequencies on the dispersive material. However, for other materials with more dispersion, the effect can be much more prominent.

# plot the permittivity at a few frequencies
freqs_plot = (freq_range[0], freq_range[1])
fig, axes = plt.subplots(1, len(freqs_plot), tight_layout=True, figsize=(12, 4))
for ax, freq_plot in zip(axes, freqs_plot):
    sim.viz_eps_2D(normal='x', cbar=True, frequency=freq_plot, ax=ax)

We can also take a look at the source to make sure it’s defined correctly over our frequency range of interst.

# Check probe and source
ax1, ax2 = sim.viz_source(source)
ax1.set_xlim(0, 1e-13)
# ax2.plot(freq_range, [1, 1])
ax2.fill_between(freq_range, [0,0], [1, 1], alpha=0.1, color='r')
ax2.legend(('source', 'measure'))

Run the simulation

We will submit the simulation to run as a new project.

# Submit a project to the cluster
project = web.new_project(sim.export(), task_name='docs_dispersion')
Uploading the json file...
Project 'docs_dispersion' status: success...

Postprocess and Plot

Once the simulation has completed, we can download the results and load them into the simulation object.

print('downloading data...')
web.download_results(project['taskId'], target_folder='output')
downloading data...
Applying source normalization to all frequency monitors using source index 0.

Now, we compute the transmitted flux and plot the transmission spectrum.

# Retrieve the power flux through the monitor plane.
transmission = sim.data(monitor)['flux']
plt.plot(monitor_lambdas, transmission, color='k')
plt.xlabel('wavelength (um)')
plt.ylabel('transmitted flux')

To get power transmission, we need to do a normalizing run without any slab and divide by that result

sim_norm = td.Simulation(
    center = (0, 0, 0),
    size = sim_size,
    mesh_step = 1/res,
    structures = [],
    sources = [source],
    monitors = [monitor],
    run_time = t_stop,
    pml_layers = (0, 0, 15)
project = web.new_project(sim_norm.export(), task_name='docs_dispersion_norm')
print('downloading results...')
web.download_results(project['taskId'], target_folder='output')
transmission_norm = sim_norm.data(monitor)['flux']
Initializing simulation...
Mesh step (micron): [6.67e-03, 6.67e-03, 6.67e-03].
Simulation domain in number of grid points: [150, 150, 1140].
Total number of computational grid points: 2.56e+07.
Total number of time steps: 21651.
Estimated data size (GB) of monitor monitor: 0.0000.
Uploading the json file...
Project 'docs_dispersion_norm' status: success...

downloading results...
Applying source normalization to all frequency monitors using source index 0.
plt.plot(monitor_lambdas, transmission, label='with structure')
plt.plot(monitor_lambdas, transmission_norm, label='no structure')
plt.plot(monitor_lambdas, transmission / transmission_norm, 'k--', label='normalized')
plt.xlabel('wavelength (um)')
plt.ylabel('fraction of transmitted power (normalized)')

We see that since the flux monitor already takes the source power into account, the normalizing run doens’t affect the calculation much, but it is prudent to include for broadband calculations.

Analytical Comparison

We will use a transfer matrix method (TMM) code to compare tidy3d transmission to a semi-analytical result.

# install and import TMM package
!pip install -q tmm
import tmm
# prepare list of thicknesses including air boundaries
d_list = [np.inf] + t_slabs + [np.inf]

# convert the complex permittivities at each frequency to refractive indices
n_list1 = np.sqrt(mat1.epsilon(monitor_freqs))
n_list2 = np.sqrt(mat2.epsilon(monitor_freqs))
n_list3 = np.sqrt(mat3.epsilon(monitor_freqs))
n_list4 = np.sqrt(mat4.epsilon(monitor_freqs))

# 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, n_list1[i], n_list2[i], n_list3[i], n_list4[i], 1]

    # get transmission at normal incidence
    T = tmm.coh_tmm('s', n_list, d_list, 0, lam)['T']
plt.plot(monitor_lambdas, transmission_tmm, label='TMM')
plt.plot(monitor_lambdas, transmission / transmission_norm, 'k--', label='Tidy3D')
plt.xlabel('wavelength ($\mu m$)')