In this lecture, we show how to use FDTD to solve a basic EM problem involving the transmission of light though a slab of material.

- How to set up the simulation using periodic boundary conditions and a plane wave source.

- How to use a broadband pulse to compute the normalized transmission over a range of wavelengths.

- How the FDTD results agree with those computed with analytical methods.

Additional information:
This Lecture was updated in Apr 14, 2022

This tutorial will walk you through computing transmission spectrum of a dielectric slab using Tidy3D.

For more details, refer to the Tidy3d documentation.

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

First, let us define some basic parameters setting up the slab and other simulation properties.

In [2]:

```
# Wavelength and frequency range
freq_range = (200e12, 400e12)
lambda_range = (td.constants.C_0/freq_range[1], td.constants.C_0/freq_range[0])
freq0 = np.sum(freq_range)/2
Nfreq = 301
# 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
lambda0 = td.C_0 / freq0
bandwidth = 0.18 # bandwidth of source in units of delta frequency. 0.38 for broadband
freqw = bandwidth * (freq_range[1] - freq_range[0])
t_stop = 100 / freq0
# Thickness and refractive index of slab
t_slab = 0.5
n_slab = 3.5
mat_slab = td.Medium(permittivity=n_slab**2, name='silicon')
# Grid resolution (cells per um)
# dl = lambda_range[0] / 30 / n_slab
dl = 6e-3
# space between slabs and sources and PML
spacing = 1 * lambda_range[-1]
# simulation size
sim_size = Lx, Ly, Lz = (2.0, 2.0, 2*spacing + t_slab)
```

First, let's use the analytical expresion for slab transmission to see what we should expect.

In [3]:

```
def slab_analytical(d, n, wvl):
""" computes transmision as a function of slab thickness (d), refractive index (n), and wavelength (wvl). """
rho = (n-1)/(n+1)
t = ((1+rho)*(1-rho)*np.exp(-2j*np.pi*n*d/wvl)) / (1 - rho**2*np.exp(-4j*np.pi*n*d/wvl))
return np.abs(t)**2
transmission_analytical = slab_analytical(t_slab, n_slab, monitor_lambdas)
```

In [4]:

```
plt.figure()
plt.plot(monitor_freqs / 1e12, transmission_analytical, 'k', label='analytical')
plt.xlabel('frequency (THz)')
plt.ylabel('Transmitted')
plt.legend()
plt.show()
```

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

First, we define the multilayer stack structure.

In [5]:

```
slab = td.Structure(
geometry=td.Box(
center=(0, 0, -Lz/2 + spacing + t_slab/2),
size=(td.inf, td.inf, t_slab),
),
medium=mat_slab,
name='slab',
)
```

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.

In [6]:

```
# 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 just past (towards positive z) of the stack.

In [7]:

```
# 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',
)
```

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.uniform(dl=dl),
structures = [slab],
sources = [source],
monitors = [monitor],
run_time = t_stop,
boundary_spec = td.BoundarySpec.pml(z=True),
normalize_index = None,
)
```

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.

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

In [10]:

```
# plot the permittivity
freqs_plot = (freq_range[0], freq_range[1])
fig, ax = plt.subplots(1, tight_layout=True, figsize=(8, 4))
sim.plot_eps(x=0, freq=None, ax=ax)
plt.show()
```

In [11]:

```
times = np.arange(0, sim.run_time, sim.dt)
amps = np.real(sim.sources[0].source_time.amp_time(times))
```

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

In [12]:

```
# Check probe and source
f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(8, 4))
ax1 = sim.sources[0].source_time.plot(times=np.linspace(0, sim.run_time, 1001), ax=ax1)
ax1.set_xlim(0, sim.run_time)
ax1.legend(('source amplitude',))
ax2 = sim.sources[0].source_time.plot_spectrum(times=np.linspace(0, sim.run_time, 1001), val='abs', ax=ax2)
fill_max = 30e-16 # 10e-16
ymax = 50e-16 # 20e-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()
```

To compare the results with and without the slab, let's make a copy of our simulation and remove the slab. This will be our reference for normalization.

In [13]:

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

We will submit both simulations to run.

In [14]:

```
sim_data0 = web.run(sim0, task_name='lecture02_slab_narrowband_normalization', path=f'data/data0_narrowband.hdf5', verbose=True)
sim_data = web.run(sim, task_name='lecture02_slab_narrowband_transmission', path=f'data/data_narrowband.hdf5', verbose=True)
```

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.

In [15]:

```
# 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(2, 1, figsize=(4, 7), tight_layout=True)
transmission0.plot(ax=ax1, label='no slab')
transmission.plot(ax=ax1, label='with slab')
transmission_normalized.plot(ax=ax2)
ax1.legend()
ax1.set_title('raw transmission')
ax2.set_title('normalized transmission')
plt.show()
```

We notice that at the edges of the frequency spectrum, the results look a bit strange.

In [16]:

```
plt.figure()
plt.plot(monitor_freqs/1e12, transmission_analytical, 'k', label='analytical')
plt.plot(monitor_freqs/1e12, transmission_normalized, 'b--', label='FDTD')
plt.xlabel('frequency (THz)')
plt.xlim([200, 400])
plt.ylabel('Transmitted')
plt.legend()
plt.show()
```

And also dont match the analytical.

Those results show some strange behavior at the edge frequencies. As we discuss in the lecture, it is likely due to the source spectrum not having enough power in those regions. So let's increase the source bandwidth and try again.

In [17]:

```
bandwidth = 0.38
freqw = bandwidth * (freq_range[1] - freq_range[0])
# make a copy of the original simulation, but change the frequency width of the source
bw_source = source.copy(update={
'source_time':source.source_time.copy(update={'fwidth':freqw})
}
)
sim_bw = sim.copy(update={'sources':[bw_source]})
# make the normalization simulation for this
sim_bw0 = sim_bw.copy(update={'structures':[]})
```

In [18]:

```
# Check probe and source, looks now like the source spectrum fully covers the green square
f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(8, 4))
ax1 = sim_bw.sources[0].source_time.plot(times=np.linspace(0, sim_bw.run_time, 1001), ax=ax1)
ax1.set_xlim(0, sim_bw.run_time)
ax1.legend(('source amplitude',))
ax2 = sim_bw.sources[0].source_time.plot_spectrum(times=np.linspace(0, sim_bw.run_time, 1001), val='abs', ax=ax2)
fill_max = 15e-16
ymax = 20e-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()
```

In [19]:

```
sim_bw_data0 = web.run(sim_bw0, task_name='lecture02_slab_broadband_normalizatio', path=f'data/data0_broadband.hdf5')
sim_bw_data = web.run(sim_bw, task_name='lecture02_slab_broadband_transmission', path=f'data/data_broadband.hdf5')
```

In [20]:

```
# Retrieve the power flux through the monitor plane.
transmission0 = sim_bw_data0['flux'].flux
transmission = sim_bw_data['flux'].flux
transmission_normalized = transmission / transmission0
f, (ax1, ax2) = plt.subplots(2, 1, figsize=(4, 7), tight_layout=True)
transmission0.plot(ax=ax1, label='no slab')
transmission.plot(ax=ax1, label='with slab')
transmission_normalized.plot(ax=ax2)
ax1.legend()
ax1.set_title('raw transmission')
ax2.set_title('normalized transmission')
plt.show()
```

In [21]:

```
plt.figure()
plt.plot(monitor_freqs/1e12, transmission_analytical, 'k', label='analytical')
plt.plot(monitor_freqs/1e12, transmission_normalized, 'b--', label='FDTD')
plt.xlabel('frequency (THz)')
plt.xlim([200, 400])
plt.ylabel('Transmitted')
plt.legend()
plt.show()
```

These results look much better.

In [ ]:

```
```

FDTD 101: Lecture 2

This is a second video on introduction to the finite difference time domain method or the FDTD method. I'm Shanhui Fan from Flexcompute.

In the last video we gave an introduction to the basic idea of the FDTD method, and showed that you can generate a very nice movie visualizing how electromagnetic waves propagate inside a vacuum region.

In this video, we're going to take a step further and to show how you can use the FDTD method to get information that you usually care about in device design. And here we're showing a particularly simple example where we're going to try to get a transmission spectrum through a dielectric slab. In our example, the dielectric slab is made of silicon with an index of 3.5 in the infrared wavelength range and with a thickness of half a micron. We will be sending light propagating normally instant upon the silicon slab along the Z axis, and try to compute the transmission. For this simple system the transmission spectrum can be computed analytically and shown here on the right is the analytic result. What we’d like to show you is how to reproduce this analytic result with the finite difference time domain method and in doing so illustrate some of the additional thinking and perhaps subtleties about the setup of the FDTD method.

To start any FDTD simulation, the first thing that you need to do is to set up the computational cell or the computational domain. In our case, we will put the silicon slab in the middle of the computational domain. We will add some free space below and above the silicon slab so that we can accommodate the source and the monitor on either side of the silicon slab. We choose 1.5 micron which corresponds to about one and a half wavelengths of free space below and above the silicon slab. In the XY direction, we choose the size of the computational domain to be 2 microns. For this problem, it is an overkill. In fact, because the silicon slab is uniform in the x and y dimension, you can in fact, choose a much smaller computational domain size, maybe even 0.1 or 0.2 micron would have worked. However, in the next video we're going to use the same computational domain to simulate the photonic crystal slab structure, which has more complexity in the in-plane direction. So the xy dimension here is chosen. We choose a discretization so that each grid corresponds to a linear dimension of about 10 nanometer and that turns out to be sufficient to get to the accuracy that we need. As I mentioned in my last video for any finite difference time domain simulation, you will need to choose the boundary condition that truncates the computational domain.

In our case, we put perfectly matched layers along the z direction and these are artificial absorbers used to absorb incident plane waves to simulate an open structure along the z direction. In the xy dimensions as I mentioned the slab is uniform. And therefore, we put in periodic boundary conditions as indicated here. These periodic boundary conditions are useful for simulating structures with infinite extent interacting within the incident plane wave.

Now let me comment on the source. Source is used to excite the electromagnetic field inside the computational domain. In the FDTD simulation, it consists of distribution of oscillating dipoles on a plane in our case. And we will have all the dipoles to have the same magnitude and they oscillate in phase in order to set up an incident plane wave.

One of the very important capabilities in FDTD simulation is that you can set up a pulsed source so that you can get a broad-band response of a given structure, in other words, the response of a given structure over a wide range of frequencies in a single simulation. In our case, we set up a source amplitude in time so that the current actually oscillates at a carrier frequency, but has a gaussian envelope in time. And this is the corresponding spectrum of the source, you can see that with this time dependent source the spectrum covers a frequency range of about a hundred terahertz around the 300 terahertz carrier frequency. We are using a power source to generate a broad-band input in order to determine the response of the structure over a broad bandwidth.

To determine the transmission, we look at the electromagnetic fields on the monitor plane on the other side of the slab, and we compute the Poynting vector flux that passes through the monitor plane as a function of frequency. This allows us to compute for a given source, how much power is actually transmitted on the other side of the structure.

To convert this into a transmission spectrum, however, we need to figure out the amount of power or intensity that's in the incident wave. Therefore, in many of these simulations we typically do two calculations. In the first calculation, we would use the same computational domain with the source and the monitor, but without the slab. We're simulating a vacuum and we look at how the Poynting vector flux looks on the monitor plane as a functional frequency and that gives us the blue curve. Then in the second simulation, we repeat the calculation by putting it into the slab. The Poynting vector flux spectrum gives us the red curve. Once we have these two curves we can divide one against the other to get the transmission spectrum and that's shown on the bottom here. You can see that the transmission spectrum in the far-infrared looks very nice. Especially the region near the 300 terahertz frequency range is very smooth. But near 200 terahertz or 400 terahertz, you can see something strange that's going on here. By analyzing the incident wave spectrum as indicated by the blue curve, you can see that at 200 terahertz or 400 terahertz, you’re really at the wing of this incidence spectrum where the power in the incident wave is very small. In this case then, in fact, you will not get very reliable data about the transmission simply because there's not enough power in the incident wave at those frequencies. To fix this what you would usually do is to choose a pulse that's narrow enough in time. That spectrum can cover the entire frequency range that you are interested in.

In our case if you narrow the temporal duration of the pulse, you can get a source spectrum as indicated by the blue curve here, which now covers a broader range of frequencies all the way from 200 terahertz to 400 terahertz. And again, repeating the calculation with the slab, and take the ratio of these two, you get a much nicer transmission spectrum as indicated here. And you can see that the curve in this case is a lot smoother especially at the wing of the instance spectrum for example from 200 terahertz to 400 terahertz.

You can compare this result with the analytic results. In this plot, the red curve is from FDTD and the dash line is from analytic results. You can see essentially perfect agreement between FDTD and analytic results. This of course is a relatively simple FDTD calculation. But if you are interested in learning FDTD, I would strongly encourage you to actually go through the calculation with the script that we provide. I think you will learn a lot on how you actually go about setting up the FDTD simulation and some of the thinking that's behind it in order to get reliable results.