TIDY3D
LEARNING CENTER

Lectures

Tutorials

Applications

Applying FDTD to Photonic Crystal Slab Simulation

By Weiliang Jin, Zongfu Yu and Shanhui Fan

In this lecture, we show how to apply FDTD to solve a slightly more complicated EM problem involving a photonic crystal slab that supports guided resonance.

- Explain the importance of setting up sufficient runtime in computing transmission spectrum, so that the field of the long-lived mode can decay to a negligible value.
- Show the long decay tail of the mode in the time domain, the transmission spectrum that exhibits Fano-resonant lineshape.

Share On:

Visit public projects:

Download .ipynb View Presentation Slides

Additional information: This was updated in Apr 13, 2022

Tutorial 3: applying FDTD to photonic crystal slab simulation

In this tutorial, we show how to apply FDTD to solve a slightly more complicated EM problem involving a photonic crystal slab that supports guided resonance. In particular, we reproduce the findings of Fan et al. (2003), which is linked here.

For more details, refer to the Tidy3d documentation.

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

td.config.local_cache.enabled = False

Next, let us define some key parameters.

runtime = 100.0  # in unit of 1/frequency bandwidth of the source
dPML = 1.0  # space between PhC slabs and PML, in unit of longest wavelength of interest

Let us define some basic parameters.

# Wavelength and frequency range
freq_range = (75e12, 135e12)
lambda_range = (td.constants.C_0 / freq_range[1], td.constants.C_0 / freq_range[0])
freq0 = np.sum(freq_range) / 2

# central frequency, frequency pulse width
lambda0 = td.C_0 / freq0
width = 0.3
freqw = width * (freq_range[1] - freq_range[0])

# runtime
t_stop = runtime / freqw
print(f"Total runtime <= {t_stop * 1e12} ps")

# frequencies and wavelengths of monitor
Nfreq = 1001
monitor_freqs = np.linspace(freq_range[0], freq_range[1], Nfreq)
monitor_lambdas = td.constants.C_0 / monitor_freqs

# PhC slab
period = 1.0  # um, period along x- and y- direction
t_slab = 0.55  # um, thickness
ep_slab = 12  # permittivity
mat_slab = td.Medium(permittivity=ep_slab, name="silicon")

# PhC air hole
radius = 0.2  # um, radius
ep_hole = 1  # permittivity
mat_hole = td.Medium(permittivity=ep_hole, name="air")

# Grid size # um
dl = (
    lambda_range[0] / 30 / np.sqrt(ep_slab)
)  # 30 grids per smallest wavelength in medium
print(f"dl = {dl * 1000} nm")

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

# simulation size
sim_size = Lx, Ly, Lz = (period, period, 2 * spacing + t_slab)

Total runtime <= 5.555555555555555 ps
dl = 21.368550163866615 nm

Create Simulation

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

First, we define the unit cell of PhC, made of a high-index slab and an air hole.

slab = td.Structure(
    geometry=td.Box(
        center=(0, 0, 0),
        size=(td.inf, td.inf, t_slab),
    ),
    medium=mat_slab,
    name="slab",
)

hole = td.Structure(
    geometry=td.Cylinder(
        center=(0, 0, 0),
        axis=2,
        radius=radius,
        length=t_slab,
    ),
    medium=mat_hole,
    name="air hole",
)

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=freq0, fwidth=freqw),
    size=(td.inf, td.inf, 0),
    center=(0, 0, -Lz / 2 + spacing / 2),
    direction="+",
    pol_angle=0,
    name="planewave",
)

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

# 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 / 2),
    size=(td.inf, td.inf, 0),
    freqs=monitor_freqs,
    name="flux",
)

Additionally, we also place a time monitor in the middle of the field monitor plane to monitor the time evolution of the field.

# We are interested in measuring the time evolution at a single point in the field monitor plane
monitor_time = td.FieldTimeMonitor(
    center=(0, 0, Lz / 2 - spacing / 2),
    size=(0, 0, 0),
    fields=["Ex", "Ey", "Ez", "Hx", "Hy", "Hz"],
    interval=1,
    name="time",
)

Now it is time to define the simulation object.

sim = td.Simulation(
    center=(0, 0, 0),
    size=sim_size,
    grid_spec=td.GridSpec.uniform(dl=dl),
    structures=[slab, hole],
    sources=[source],
    monitors=[monitor, monitor_time],
    run_time=t_stop,
    shutoff=1e-7,
    boundary_spec=td.BoundarySpec.pml(z=True),
    normalize_index=None,
)

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, 2, tight_layout=True, figsize=(10, 4))
sim.plot(z=0.0, ax=ax[0])
sim.plot(x=0, freq=freq0, ax=ax[1])
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
f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(8, 4))
plot_time = 0.2e-12
ax1 = sim.sources[0].source_time.plot(times=np.linspace(0, plot_time, 1001), ax=ax1)
ax1.set_xlim(0, plot_time)
ax1.legend(("source amplitude",))
ax2 = sim.sources[0].source_time.plot_spectrum(
    times=np.linspace(0, sim.run_time, 10001), val="abs", ax=ax2
)
fill_max = 45e-16
ymax = 60e-16
ax2.fill_between(
    freq_range, [-0e-16, -0e-16], [fill_max, fill_max], alpha=0.4, color="g"
)
ax2.legend(("source spectrum", "measurement"))
ax2.set_ylim(-1e-16, ymax)
plt.show()

Finally, we define a simulation object without PhC slab to normalize transmission.

sim0 = sim.copy(update={"structures": []})

Running Simulation

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

batch = {"lecture03_PhC_normalization": sim0, "lecture03_PhC_transmission": sim}

batch_data = web.run_async(simulations=batch, path_dir="data")

sim_data0 = batch_data["lecture03_PhC_normalization"]
sim_data = batch_data["lecture03_PhC_transmission"]

16:04:21 UTC Started working on Batch containing 2 tasks.

             Maximum FlexCredit cost: 0.081 for the whole batch.

             Use 'Batch.real_cost()' to get the billed FlexCredit cost after    
             completion.

16:04:23 UTC Batch complete.

Post Run Analysis

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.

transmission0 = sim_data0["flux"].flux
transmission = sim_data["flux"].flux
transmission_normalized = transmission / transmission0

f, (ax1, ax2) = plt.subplots(1, 2, figsize=(7, 4), tight_layout=True)
transmission0.plot(ax=ax1, label="no PhC")
transmission.plot(ax=ax1, label="with PhC")
ax1.set_xlim((freq_range[0], freq_range[1]))
transmission_normalized.plot(ax=ax2)
ax2.hlines(1.0, freq_range[0], freq_range[1], linestyles="dashed", color="black")
ax2.hlines(0.0, freq_range[0], freq_range[1], linestyles="dashed", color="black")
ax2.set_xlim((freq_range[0], freq_range[1]))
ax1.legend()
ax1.set_title("raw transmission")
ax2.set_title("normalized transmission")
plt.show()

To compare to the results in Fig.2 of the paper, let’s plot the transmission as a function of frequency in the unit of c/period

plt.figure()
plt.plot(
    monitor_freqs * period / td.C_0,
    transmission_normalized,
    label="normalized transmission",
)
plt.xlabel("frequency (c/period)")
plt.xlim([0.25, 0.45])
plt.hlines(1.0, 0.25, 0.45, linestyles="dashed", color="black")
plt.hlines(0.0, 0.25, 0.45, linestyles="dashed", color="black")
plt.ylabel("Transmitted")
plt.show()

Finally, we look into the time evolution of Ex at the center of the flux monitor.

time_data = sim_data["time"]
time_data0 = sim_data0["time"]
fig, ax = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True)
ylim = 2
xlim = t_stop
# with PhC
Ex = time_data.Ex
Ex.plot(ax=ax[1])
ax[1].set_ylim([-ylim, ylim])
ax[1].set_xlim([0, xlim])
ax[1].set_ylabel("$E_x(t)$ [V/m]")
ax[1].set_title("with PhC")
# without PhC
Ex0 = time_data0.Ex
Ex0.plot(ax=ax[0])
ax[0].set_ylabel("$E_x(t)$ [V/m]")
ax[0].set_xlim([0, xlim])
ax[0].set_ylim([-ylim, ylim])
ax[0].set_title("without PhC")

plt.show()

Without PhC, the field decays rapidly; with PhC, a long decay tail is clearly visible.

Download Photonic Crystal Slab Simulation FDTD

TIDY3D LEARNING CENTER