TIDY3D
LEARNING CENTER

Lectures
Tutorials
Applications

Example Library / Mediums / Modeling dispersive materials

Modeling dispersive materials

Sign up now to simulate .ipynb

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.

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

First, we define basic simulation-related parameters. To simplify the initialization and management of frequency-related variables, we use the convenience class FreqRange.

# Wavelength interval and number of sampled frequencies
wvl_min = 0.5
wvl_max = 1.5
Nfreq = 333

# initialize frequency range, wavelengths of monitor
freq_range = td.FreqRange.from_wvl_interval(wvl_min=wvl_min, wvl_max=wvl_max)
monitor_lambdas = freq_range.wvls(num_points=Nfreq, spacing="uniform_freq")[::-1]

# central frequency, frequency pulse width and total running time
t_stop = 100 / freq_range.freq0

# Thicknesses of slabs
t_slabs = [0.5, 0.2, 0.4, 0.3]  # um

# Grid resolution (min steps per wavelength in a material)
res = 40

# space between slabs and sources and PML
spacing = wvl_max

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

  1. Simple, lossy dielectric defined by a real-valued relative permittivity, and DC conductivity.
  2. 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.
  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. 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
  4. Dispersive material imported from our predefined 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.

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

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

# Here we define the planewave source, placed just in advance (towards negative z) of the slab
source = td.PlaneWave(
    source_time=td.GaussianPulse(freq0=freq_range.freq0, fwidth=0.3 * freq_range.fwidth),
    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.

# 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=freq_range.freqs(num_points=Nfreq, spacing="uniform_freq"),
    name="flux",
)

Next, define the boundary conditions to use PMLs along z and the default periodic boundaries along x and y

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.

sim = td.Simulation(
    center=(0, 0, 0),
    size=sim_size,
    grid_spec=td.GridSpec.auto(min_steps_per_wvl=res),
    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.

sim.plot(x=0)
plt.show()

Next, we use Simulation.plot_eps() to visualize 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 frequencies 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.freqs(num_points=2)
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 interest.

# 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.freqs(num_points=2),
    [-8e-16, -8e-16],
    [8e-16, 8e-16],
    alpha=0.4,
    color="g",
    label="measure",
)
ax2.legend()
plt.show()

Run the simulation

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

sim_data = web.run(sim, task_name="dispersion", path="data/sim_data.hdf5", verbose=True)
08:07:17 UTC Created task 'dispersion' with resource_id                         
             'fdve-23bcf743-ccf4-42f1-a950-6c3948f1c1ae' and task_type 'FDTD'.  

             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-23bcf743-ccf
             4-42f1-a950-6c3948f1c1ae'.                                         

             Task folder: 'default'.                                            







08:07:49 UTC Estimated FlexCredit cost: 0.104. 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.                                                             


















08:08:04 UTC starting up solver                                                 







             running solver                                                     












08:08:18 UTC early shutoff detected at 28%, exiting.                            













08:08:19 UTC status = postprocess                                               












             status = success                                                   













08:08:21 UTC View simulation result at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-23bcf743-ccf
             4-42f1-a950-6c3948f1c1ae'.                                         


















08:08:22 UTC Loading results from data/sim_data.hdf5                            











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.






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






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







             Created task 'docs_dispersion_norm' with resource_id               
             'fdve-d4276c7a-05d4-4325-9815-42f7c2158893' and task_type 'FDTD'.  







             View task using web UI at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-d4276c7a-05d
             4-4325-9815-42f7c2158893'.                                         







             Task folder: 'default'.                                            


















08:08:23 UTC Estimated FlexCredit cost: 0.025. 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.                     







08:08:24 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.                                                             












08:08:36 UTC status = preprocess                                                







08:08:39 UTC status = queued                                                    







08:08:57 UTC status = preprocess                                                













08:09:01 UTC starting up solver                                                 







08:09:02 UTC running solver                                                     












             early shutoff detected at 17%, exiting.                            













08:09:03 UTC status = success                                                   













             View simulation result at                                          
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-d4276c7a-05d
             4-4325-9815-42f7c2158893'.                                         


















             Loading results from data/sim_data.hdf5                            











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.






# import TMM package
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.eps_model(freq_range.freqs(num_points=Nfreq)))
n_list2 = np.sqrt(mat2.eps_model(freq_range.freqs(num_points=Nfreq)))
n_list3 = np.sqrt(mat3.eps_model(freq_range.freqs(num_points=Nfreq)))
n_list4 = np.sqrt(mat4.eps_model(freq_range.freqs(num_points=Nfreq)))

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











plt.figure()
plt.plot(monitor_lambdas, transmission_tmm, label="TMM")
plt.plot(monitor_lambdas, transmission / transmission_norm, "k--", label="Tidy3D")
plt.xlabel(r"wavelength ($\mu m$)")
plt.ylabel("Transmitted")
plt.legend()
plt.show()



























          

          

            Example library
          

          

          
          

          

            
            
              
              Previous Tutorial
            
            
            
            
              Next Tutorial