Modeling dispersive materials#
Introduction / Setup#
Here we show how to model dispersive materials in Tidy3D with an example showing transmission spectrum of a multilayer stack of slabs.
If you are new to the finite-difference time-domain (FDTD) method, we highly recommend going through our FDTD101 tutorials. For simulation examples, please visit our examples page. FDTD simulations can diverge due to various reasons. If you run into any simulation divergence issues, please follow the steps outlined in our troubleshooting guide to resolve it.
[1]:
# standard python imports
import numpy as np
import matplotlib.pyplot as plt
import tidy3d as td
from tidy3d import web
First, let us define some basic parameters.
[2]:
# 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)#
Here, we will illustrate defining materials in four different ways:
Simple, lossy dielectric defined by a real-valued relative permittivity, and DC conductivity.
Active material defined by real and imaginary part of the refractive index (\(n\)) and (\(k\)) at a given frequency. Values are exact only at that frequency, so this approach is only good for narrow-band simulations.
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. The dispersion must be negative. This is a convenient approach to incorporate weakly dispersive materials in your simulations, as the values can be taken directly from refractiveindex.info
Dispersive material imported from our pre-defined library of materials.
More complicated dispersive materials can also be defined through dispersive models like Lorentz, Sellmeier, Debye, or Drude, if the model parameters are known. Finally, arbitrary dispersion data can also be fit, which is a the subject of this tutorial.
[3]:
# simple, lossy material
mat1 = td.Medium(permittivity=4.0, conductivity=0.005)
# active material with n & k values at a specified frequency or wavelength
# note: negative k value corresponds to a gain medium; it is only allowed
# when `allow_gain` is set to be True
mat2 = td.Medium.from_nk(n=3.0, k=-0.1, freq=freq0, allow_gain=True)
# weakly dispersive material defined by dn_dwvl at a given frequency
mat3 = td.Sellmeier.from_dispersion(n=2.0, dn_dwvl=-0.1, freq=freq0)
# dispersive material from tidy3d library
mat4 = td.material_library["BK7"]["Zemax"]
# 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.
[4]:
slabs = []
slab_position = -Lz / 2 + 2 * spacing
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,
)
slabs.append(slab)
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.
[5]:
# 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),
direction="+",
pol_angle=0,
)
Here we define the field monitor, placed just past (towards positive z) of the stack.
[6]:
# 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="flux",
)
Next, define the boundary conditions to use PMLs along z and the default periodic boundaries along x and y
[7]:
boundary_spec = td.BoundarySpec(
x=td.Boundary.periodic(), y=td.Boundary.periodic(), z=td.Boundary.pml()
)
Now it is time to define the simulation object.
[8]:
sim = td.Simulation(
center=(0, 0, 0),
size=sim_size,
grid_spec=td.GridSpec.auto(min_steps_per_wvl=40),
structures=slabs,
sources=[source],
monitors=[monitor],
run_time=t_stop,
boundary_spec=boundary_spec,
)
Plot The Structure#
Let’s now plot the permittivity profile to confirm that the structure was defined correctly.
First 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.
[9]:
sim.plot(x=0)
plt.show()
Next, we use Simulation.plot_eps() to vizualize the permittivity of the stack. However, because the stack contains dispersive materials, we need to specify the freq 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.
[10]:
# 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.plot_eps(x=0, freq=freq_plot, ax=ax)
plt.show()
We can also take a look at the source to make sure it’s defined correctly over our frequency range of interst.
[11]:
# Check probe and source
ax1 = sim.sources[0].source_time.plot(times=np.linspace(0, sim.run_time, 1001))
ax1.set_xlim(0, 1e-13)
ax2 = sim.sources[0].source_time.plot_spectrum(times=np.linspace(0, sim.run_time, 1001))
ax2.fill_between(
freq_range, [-8e-16, -8e-16], [8e-16, 8e-16], alpha=0.4, color="g", label="mesure"
)
ax2.legend()
plt.show()
Run the simulation#
We will submit the simulation to run as a new project.
[12]:
sim_data = web.run(sim, task_name="dispersion", path="data/sim_data.hdf5", verbose=True)
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 856a276d-7d85-4f17-b50c-be3e3f57777bv1'.
Postprocess and Plot#
Once the simulation has completed, we can download the results and load them into the simulation object.
Now, we compute the transmitted flux and plot the transmission spectrum.
[13]:
# Retrieve the power flux through the monitor plane.
transmission = sim_data["flux"].flux
plt.plot(monitor_lambdas, transmission, color="k")
plt.xlabel("wavelength (um)")
plt.ylabel("transmitted flux")
plt.show()
In Tidy3D, results are normalized by default. In some cases, and largely depending on the required accuracy, a normalizing run may still be needed. Here, we show how to do such a normalizing run by simulating an empty simulation with the exact same source and monitor but none of the structures.
[14]:
sim_norm = sim.copy(update={"structures": []})
sim_data_norm = web.run(
sim_norm, task_name="docs_dispersion_norm", path="data/sim_data.hdf5", verbose=True,
)
transmission_norm = sim_data_norm["flux"].flux
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 66335a77-705d-42a5-be83-4aac692b3f9fv1'.
[15]:
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.legend()
plt.xlabel("wavelength (um)")
plt.ylabel("fraction of transmitted power (normalized)")
plt.show()
We see that since the flux monitor already takes the source power into account, the normalizing run has no visible effect on the results.
Analytical Comparison#
We will use a transfer matrix method (TMM) code to compare Tidy3D’s simulated transmission to a semi-analytical result.
[16]:
# import TMM package
import tmm
[17]:
# 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.eps_model(monitor_freqs))
n_list2 = np.sqrt(mat2.eps_model(monitor_freqs))
n_list3 = np.sqrt(mat3.eps_model(monitor_freqs))
n_list4 = np.sqrt(mat4.eps_model(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"]
transmission_tmm.append(T)
[18]:
plt.figure()
plt.plot(monitor_lambdas, transmission_tmm, label="TMM")
plt.plot(monitor_lambdas, transmission / transmission_norm, "k--", label="Tidy3D")
plt.xlabel("wavelength ($\mu m$)")
plt.ylabel("Transmitted")
plt.legend()
plt.show()
[ ]: