Perfectly matched layer (PML) is commonly used to truncate unbounded computational region, since an ideal PML can completely absorb the incoming waves from all angles of the incidence without any reflection. In this lecture, we explain the basic idea behind PML, and show how to characterize the performance of PML in FDTD simulations.

- Explain the concept of PML through coordinate transformation.

- Show that the performance of PML should be characterized by the amplitude rather than the intensity of parasitic reflection, and present a simple approach to compute parasitic reflection at many frequencies.

- Illustrate typical signature in the spectrum when the parasitic reflection at the boundary is significant.

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

In this notebook, we show how to measure the parasitic reflection from the perfectly matched layer (PML) boundary in order to benchmark the performance of PML. As a simple illustration, we consider a 1D system where a time monitor is placed between a broadband planewave source and PML bounday. By visualizing the evolution of the field at the time monitor, one can clearly distinguish between the incident pulse passing through the time monitor and the reflected pulse. A subsequent Fourier analysis can tell the reflection amplitude at the PML boundary over a broad range of frequencies.

Finally, taking Fabry-Parot cavity as an example, we illustrate the impact of reflection from PML.

In [1]:

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

Next, let us define some key parameters.

In [2]:

```
# The setup of the 1D simulation:8
# | | | | | | | |
# |`PML`| dPML |`source`| distance |`monitor`| distance |`PML`|
# | | | | | | | |
# Number of PML layers
Npml = 2 # the default number of layers is 12
# sufficiently large distance between the source and the time monitor
distance = 100 #um
# distance between the source and PML
dPML = 1 #um
# Wavelength and frequency range of the source
freq0 = 300e12 #Hz, central frequency
fwidth = freq0/5 #frequency pulse width
lambda0 = td.C_0 / freq0
# resolution
resolution = 20
dl = lambda0/resolution
# runtime
t_stop = 5e-12
# simulation domain size
sim_size = Lx, Ly, Lz = (dl, dl, 2*distance+dPML)
```

In [3]:

```
# Uniform grid
grid_spec = td.GridSpec.uniform(dl=dl)
# Planewave source
src_z = dPML - Lz/2
source = td.PlaneWave(
center=(0,0,src_z),
size=(td.inf, td.inf, 0),
source_time = td.GaussianPulse(
freq0=freq0,
fwidth=fwidth,
),
direction='+',
name='Gaussian pulse',
)
# Time monitor at a single point
tmnt_pt = td.FieldTimeMonitor(
center=(0, 0, src_z + distance),
name='time',
size=(0, 0, 0),
)
# PML boundary
pml = td.Boundary.pml(num_layers=Npml)
bspec = td.BoundarySpec(x=td.Boundary.periodic(), y = td.Boundary.periodic(), z=pml)
```

Now it is time to define the simulation object.

In [4]:

```
sim = td.Simulation(
center = (0, 0, 0),
size = sim_size,
grid_spec = grid_spec,
structures = [],
sources = [source],
monitors = [tmnt_pt],
run_time = t_stop,
boundary_spec = bspec,
shutoff = 0,
)
```

In [5]:

```
sim_data = web.run(sim, task_name='lecture06_pml_reflection', path=f'data/data_pml_reflection.hdf5')
```

We visualize the time evolution of the field at the monitor plane. Two pulses are clearly visible in the figure. The amplitude of the pulse appearing first is much larger than the one appearing later. To identify the pulses, we compute analytically the time it takes for the center of the pulse to reach the monitor (red dashed line) and for it to get reflected and reach the monitor again (blue dashed line).

In [6]:

```
plt.rcParams.update({'font.size': 16})
fig, ax = plt.subplots(1,figsize=(7,4), tight_layout=True)
speed = td.C_0
time_reach_monitor = source.source_time.offset/fwidth/2/np.pi + distance/speed
time_reflection = time_reach_monitor + 2*distance/speed
ax.vlines(time_reach_monitor*1e12,-600,600,colors='r',linestyle='dashed')
ax.vlines(time_reflection*1e12,-600,600,colors='b',linestyle='dashed')
time_data = sim_data['time']
Ex = time_data.Ex
Ex_values = np.squeeze(Ex.values)
ax.plot(Ex.t*1e12,Ex_values)
ax.set_xlabel('time (ps)')
ax.set_ylabel('$E_x$')
ax.set_xlim(0,1.5)
ax.set_ylim(-600,600)
plt.show()
```

Now to compute the reflection amplitude across the frequency range spanned by the Gaussian source, let's perform Fourier analysis.

In [7]:

```
# First, let's select two time intervals of equal length. The first interval
# completely contains the original incident pulse, and the second interval
# completely contains the reflected pulse. We take the former to
# be [0,time_reach_PML), and the latter to be [time_reach_PML,2*time_reach_PML)
time_reach_PML = source.source_time.offset/fwidth/2/np.pi + 2* distance/speed
inds = np.argwhere(np.array(Ex.t)>=time_reach_PML)
Nt = int(inds[0])
```

In [8]:

```
# Now let's perform FFT in each interval, and visualize the spectrum
dt = sim.dt
df = 1/Nt/dt
fmesh = np.arange(0, Nt)*df
# original pulse
spectrum_original = 1/np.sqrt(2*np.pi)*np.fft.ifft(np.fft.ifftshift(Ex_values[:Nt]))/df
#reflected pulse
spectrum_reflection = 1/np.sqrt(2*np.pi)*np.fft.ifft(np.fft.ifftshift(Ex_values[Nt:Nt*2]))/df
# visualization of each spectrum
fig, ax = plt.subplots(1,2,figsize=(8,4), tight_layout=True)
ax[0].plot(fmesh/1e12, np.abs(spectrum_original),label='Original')
ax[0].plot(fmesh/1e12, np.abs(spectrum_reflection),label='Reflected')
ax[0].set_xlim((freq0-fwidth)/1e12,(freq0+fwidth)/1e12)
ax[0].set_ylim(0,8e-13)
ax[0].set_xlabel('Frequency (THz)')
ax[0].set_title('Spectrum')
ax[0].legend()
# visualization of reflection spectrum
ax[1].plot(fmesh/1e12, np.abs(spectrum_reflection/spectrum_original))
ax[1].set_xlim((freq0-fwidth)/1e12,(freq0+fwidth)/1e12)
ax[1].set_yscale('log')
ax[1].set_ylim(1e-6,1)
ax[1].set_title('Reflection amplitude')
ax[1].set_xlabel('Frequency (THz)')
plt.show()
```

Finally, let's examine the impact of parasitic reflection from PML boundary to a simple Fabry-Perot cavity transmission simulation illustrated in our previous Lecture 2. Most setups are similar, so please checkout the link here for more detailed explanations.

In [9]:

```
# 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
# frequencies and wavelengths of monitor
Nfreq = 301
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.4 # bandwidth of source in units of delta frequency. 0.38 for broadband
freqw = bandwidth * (freq_range[1] - freq_range[0])
t_stop = 500 / 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 = 6e-3
# space between slabs and sources and PML
spacing = 16 * lambda_range[-1]
# simulation size
sim_size = Lx, Ly, Lz = (dl, dl, 2*spacing + t_slab)
```

In [10]:

```
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)
fig, ax = plt.subplots(1,figsize=(5,4), tight_layout=True)
ax.plot(monitor_freqs / 1e12, transmission_analytical, 'k', label='analytical')
ax.set_xlabel('frequency (THz)')
ax.set_ylabel('Transmitted')
ax.legend()
plt.show()
```

In [11]:

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

In [12]:

```
# simulation with default 12 PML layers
pml = td.Boundary.pml()
bspec = td.BoundarySpec(x=td.Boundary.periodic(), y = td.Boundary.periodic(), z=pml)
sim_default = 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 = bspec,
)
# simulation with only 2 PML layers
pml_2layer = td.Boundary.pml(num_layers=2)
bspec_2layer = td.BoundarySpec(x=td.Boundary.periodic(), y = td.Boundary.periodic(), z=pml_2layer)
sim_2layer = sim_default.copy(update={'boundary_spec':bspec_2layer})
```

In [13]:

```
sim_FP_data = web.run(sim_default, task_name='lecture06_FP_default_PML', path=f'data/data_FP_default.hdf5')
sim_FP_2layer_data = web.run(sim_2layer, task_name='lecture06_FP_2PMLlayer', path=f'data/data_FP_2layer.hdf5')
```

Now, we compute the transmitted flux and plot the transmission spectrum and compare them to the analytics.

In [14]:

```
fig, ax = plt.subplots(1,figsize=(6,4), tight_layout=True)
transmission_default = sim_FP_data['flux'].flux
transmission_2layer = sim_FP_2layer_data['flux'].flux
ax.plot(monitor_freqs/1e12, transmission_analytical, 'k', label='analytical')
ax.plot(monitor_freqs/1e12, transmission_default, 'b--', label='default')
ax.plot(monitor_freqs/1e12, transmission_2layer, 'r--', label='2layer')
ax.set_xlabel('frequency (THz)')
ax.set_xlim([200, 400])
ax.set_ylim(0,1.8)
ax.set_ylabel('Transmitted')
plt.show()
```

As shown in the plot, the transmission result with the default number of PML layers agrees well with the analytics; however, when there are 2 PML layers, the reflection at the PML boundary results in lots of ripples in the transmission spectrum. Therefore, suppressing parasitic reflection from the boundary with more PML layers is extremely important for obtaining high-quality simulation result.

In [ ]:

```
```

FDTD 101: Lecture 6

Today we would continue to talk about some of the issues related to FDTD. I'm Shanhui Fan from Flexcompute.

We would talk about the issue associated with perfectly matched layer (PML) boundary condition. In FDTD we are simulating a system in a finite computational domain. So naturally, the question is, how do you truncate this domain? And in many of the simulation we would like to look at the system where the field can escape to a far field or to Infinity. At a boundary, one would need a boundary condition where there's no reflection for every single angle of instance. As an illustration of why this is important, suppose you are interested in computing the radiation pattern of a dipole in free space. Let's say you have a computational domain, you put a dipole at the center and you try to look at the field distribution as indicated by the yellow plant here. If you surround the computational domain entirely with a perfectly electric conductor boundary condition, which is a perfect mirror, then you won't have your usual dipole radiation pattern at all. But if you do the boundary condition, right, which is the perfectly matched layer (PML) boundary condition, then you get back the expected field distribution of a type of radiation pattern. The use of the boundary condition to truncate the spatial domain is extremely important.

The idea in standard FDTD simulation is to take the computational domain and then surround them by a special lossy material called perfectly matched layer (PML) to absorb incoming waves from all angles of incidence without any reflection. It is perhaps useful to go through some of the basic ideas behind this PML concept.

One way to present the PML concept is through coordinate transformation. Imagine that you have a wave propagating in vacuum along the positive X direction, if you take a snapshot at a given time, then you see a sinusoidal wave pattern extending from minus infinity to infinity. What you would like to do in PML is to put a PML layer. For example, X greater than 0 and this accomplish a coordinate transformation, such that, you transform the plane wave as was originally on the right side of the vacuum region into an exponentially decaying wave and therefore is an attenuated wave.

The detailed math for doing so can be Illustrated first by considering an 1D wave equation. We replace the usual spatial derivative operator with the derivative operator that has a stretch factor. And this stretch factor is one, and therefore you go back to the usual wave equation for x less than zero. For X less than zero you have vacuum. For X greater than 0, you put in a complex stretch factor so that the resulting wave is attenuated in the positive X direction. For this equation with a stretch factor that looks like this, this is how the solution looks like for X less than zero. You have a perfect plane wave. For X greater than 0, you have exponential decay. There's only a wave propagating along the positive X direction and there's no back reflecting wave. All the amplitude incident on the PML enters into the PML, and then gets attenuated wave. The form of the stretch factor, the sigma, is essentially a conductivity. By choosing a conductivity that is uniform, that is frequency independent. This stretch factor here is frequency dependent, and the result is that the attenuation factor here has no frequency dependency. So, in one dimension this choice would have a frequency independent attenuation. Therefore you will be able to absorb waves over a broad frequency range.

Now going from one dimension to higher dimension. Again, you can start with the wave equation and then you basically perform this coordinate transformation, along the x direction for an interface that's normal to the x direction. You can use the same stretch factor to obtain a solution that looks like this. Here is a visualization of the solution. What you see on the left in the vacuum region, is a plane wave that's incident on the interface with no reflection. And on the right, you see a wave that's attenuated inside a PML region. You can see that in this case one introduces an attenuation factor, that is only along the X Direction without touching the Y and Z coordinate. There's no attenuation along the Y and Z Direction. Therefore, a PML region is inherently anisotropic. And this anisotropy is important because it allows you to achieve zero reflection for every angle of incidence. And the choice of the attenuation only along the X direction is very important because the k_y and k_z which are conserved across the interface along interfaces, there's no attenuation. There's a matching between the vacuum region and the PML region. And this at least analytically is very attractive because you give you zero reflection and very high absorption for every angle of incidence. So here we are showing the concept of PML. The coordinate transformation concept is continuous and applied to a partial differential equation.

In FDTD, we discretize the grid. So that introduces several issues. For example, in the exact wave equation, the PML is reflectionless, even when there is a discontinuity in the conductivity profile. But reflection is in fact introduced if you solve for the discretized problem. And one of the common techniques that's being used is to taper the conductivity profile in real space. In the PML region the field decays exponentially, but you will still need a certain thickness of the PML. So that the exponential field tail reflects on the boundary of the PML doesn't strongly come back into the computational cell. You need a certain number of PML layers. There is, in fact, a very large literature formulating PML for different purposes.

Having introduced the concept of PML, I think it will be useful, maybe to talk a little bit about how these parasitic reflections from the boundaries influence the quality of your simulation. Suppose typically when you talk about these simulations, you are typically interested in, for example, intensity transmission coefficient, or intensity reflection coefficient. The desired quantity that you usually care about is the field absolute value. Now suppose you have parasitic reflection. Then the field is going to be different from the actual physical field by a perturbation as introduced by the reflection. As a result, the intensity is also going to be different from the desired intensity. And when you look at the formula here, you will see that the leading factor of the deviation from the accurate results scales as delta E, the coefficient of amplitude of the reflection coefficient. Even though you are measuring intensity, which is absolute value E square, the error in measuring the intensity due to parasitic reflection scales as the amplitude not the intensity of the reflected wave. When we characterize the reflectivity of a PML, the relevant reflectivity is the amplitude reflectivity.

Now you can numerically characterize the performance of a PML. One way to do that is you can put in a gaussian source in time to generate a wave that propagates in the computational region. Then you put a monitor plan somewhere between the gaussian source and the PML region at the end. Here is a typical plot where you show the field as a function of time. Initially you will see a pulse that passes right through the monitor region and then this pulse will keep propagating hitting the PML region and then come back. You will see the reflected pulse. Now once you have this, if you simply take those time steps corresponding to the incident pulse and those that corresponding to the reflected pulse you will be able to get the reflection coefficient.

You can now transform the input in a Fourier transformation to measure the ratio between the two then you get the reflection amplitude. So here are the simulation results for two layer or PML regions, you get a reflectivity on the order of 10%. But if you take 12 layers, it greatly improves. The reflection reduces to about 10,000 times.

One of the great things about PML is that you give angle independent attenuation. In this case here is a computation where you look at the reflection amplitude for 2, 4 and 12 layers as a function of angle of incidence, going all the way from zero to 40 degrees. You see that again the 12 layer of course has much lower reflectivity compared with the two layer case and also you can see that it is very flat as a function of angle.

So now finally I would like to directly illustrate how this kind of different reflectivity is going to influence the quality of your simulation. We'll be comparing a simulation where we use only two layers of PML compared with the case where we use 12 layers of PML. It is a simple structure: a Fabry-Perot cavity or a thin layer of silicon. And we look at the transmission coefficient for light normally instant upon the structure. You can plot the transmission as a function of frequency. The blue curve here is the result that you obtain from 12 layers of PML and this is actually almost identical to the analytic result. If you recall for a two layer PML, the reflectivity amplitude is only order of 10% or so. With that reflectivity you can see very large parasitic oscillations due to the reflection from the PML and that is indicated by the red curve here. The 12 layer, or the correct result in a way, cuts through this very large oscillation. This is a typical signature that if you have parasitic boundary reflectivity in the spectrum, you will typically see this oscillation. And that's something that if you see in your actual numerical simulation, you will have a sense of what to do. To briefly conclude, what I hope to give you a sense is the basic argument about how a PML actually works and also the fact that suppressing the parasitic reflection is extremely important for actual simulation in order to obtain high quality results.