Note: the cost of running the entire notebook is larger than 1 FlexCredit.
The direct binary search (DBS) algorithm is an iterative optimization algorithm often used in image processing and has been adapted for use in photonic design.
The steps in a DBS are typically as follows:
-
Initialization: Starting with an initial design, usually a binary representation of the photonic structure where each element (pixel or voxel) is either one material or another.
-
Evaluation: Simulating the optical performance of the current design using methods like finite-difference time-domain (FDTD) simulations.
-
Modification: Adjusting the design by flipping the state of specific elements (from one material to another) based on certain criteria, such as improving performance metrics like transmission efficiency.
-
Iteration: Repeating the evaluation and modification steps until a satisfactory design is reached or a certain number of iterations have been completed.
This process allows for the fine-tuning of photonic structures to achieve highly specific and efficient optical characteristics. The DBS algorithm is valued for its relative simplicity and effectiveness, especially in scenarios where gradient-based optimization is challenging.
This notebook demonstrates the design of an on-chip optical switch based on the phase change material Ge$_2$Sb$_2$Te$_5$ (GST) using the DBS algorithm. The design region is a square of 900 nm by 900 nm discretized into 9 by 9 pixels with a size of 100 nm. The 3 by 3 pixels at the center of the design region is GST. Other pixels can be either silicon or silica. The goal is to optimize the materials at each pixel such that the best switching performance is obtained. That is, we want to achieve the highest transmission possible when GST is in the amorphous state while the lowest transmission possible when GST is in the crystalline state. The notebook follows the work Yin, K.; Gao, Y.; Shi, H.; Zhu, S. Inverse Design and Numerical Investigations of an Ultra-Compact Integrated Optical Switch Based on Phase Change Material. Nanomaterials 2023, 13, 1643. DOI: 10.3390/nano13101643.
Tidy3D is a powerful tool for photonic design optimization due to its fast speed and high throughput. Besides DBS, we have demonstrated particle swarm optimizations of a polarization beam splitter and a bullseye cavity, genetic algorithm optimization of an on-chip reflector, and CMA-ES optimization of a waveguide S-bend. Furthermore, we also have a growing list of gradient-based adjoint optimization examples including
import numpy as np
import matplotlib.pyplot as plt
import tidy3d as td
import tidy3d.web as web
Simulation Setup¶
The simulation wavelength range is 1535 nm to 1565 nm.
lda0 = 1.55 # central wavelength
freq0 = td.C_0 / lda0 # central frequency
ldas = np.linspace(1.535, 1.565, 11) # wavelength range
freqs = td.C_0 / ldas # frequency range
fwidth = 0.5 * (np.max(freqs) - np.min(freqs)) # width of the source frequency range
For simplicity, we will use the silicon and oxide media directly from the material library. The GST has two states: the crystalline state with a higher loss and an amorphous state with a lower loss. Since our wavelength range of interest is relatively narrow, we use the from_nk method to define the material responses of GST in two states. It's important to note that the defined refractive index is strictly accurate only at the given frequency.
Si = td.material_library["cSi"]["Palik_Lossless"]
SiO2 = td.material_library["SiO2"]["Palik_Lossless"]
aGST = td.Medium.from_nk(n=3.98, k=0.0244, freq=freq0)
cGST = td.Medium.from_nk(n=6.49, k=1.054, freq=freq0)
The straight waveguide has a width of 400 nm and a thickness of 220 nm. The size of the pixels in the design region is 100 nm by 100 nm. There are 9 by 9 pixels in total but the 3 by 3 pixels at the center are fixed to GST. In the current version, we only support Nx, Ny, and N_gst being odd numbers.
w = 0.4 # width of the waveguide
t = 0.22 # thickness of the silicon
l = 1 # length of the waveguide in the simulation
Px = Py = 0.1 # pixel sizes in the x and y directions
Nx = 9 # number of pixels in the x direction
Ny = 9 # numbre of pixels in the y direction
N_gst = 3 # number of pixels for GST
buffer = 0.8 # buffer spacing
Now we want to characterize a design by a parameter vector $(x_1, x_2, ..., x_n)$, where $x_i$ represents the material at the i-th pixel and takes the value of either 0 (silica) or 1 (silicon). Due to the symmetry, we only focus on symmetric designs. Since the GST pixels are fixed, we have a total of (Ny * Nx - N_gst * N_gst) // 2 + Nx // 2 - 1 = 39 tuneable pixels. That is, the length of the parameter vector is 39 in this case.
For easier structure construction, we define a helper function parse_to_2d that takes the parameter vector and parses it into a 2D array representing the design. With that, we define another helper function create_design to define the structures in the design region as well as the input and output straight waveguides.
def parse_to_2d(array_1d, Nx, Ny, N_gst):
# ensure Ny, Nx, and N_gst are odd
if Ny % 2 == 0 or Nx % 2 == 0 or N_gst % 2 == 0:
raise ValueError("Nx, Ny, and N_gst must be odd numbers.")
# calculate the number of changeable elements in the top half (including the middle row)
num_elements_half = (Ny * Nx - N_gst * N_gst) // 2 + Nx // 2 - 1
if len(array_1d) != num_elements_half:
raise ValueError(f"1D array must have {num_elements_half} elements.")
array_2d = np.zeros((Ny, Nx))
idx_1d = 0
for i in range(Ny):
for j in range(Nx):
# skip the fixed center
if (Ny - N_gst) // 2 <= i < (Ny + N_gst) // 2 and (Nx - N_gst) // 2 <= j < (
Nx + N_gst
) // 2:
continue
# Fill up to and including the middle row
if i <= Ny // 2:
array_2d[i, j] = array_1d[idx_1d]
idx_1d += 1
# Mirror the element to the symmetrical position, but not for the middle row
if i < Ny // 2:
array_2d[Ny - 1 - i, j] = array_1d[idx_1d - 1]
return array_2d
def create_design(pixels):
# parse 1d parameter vector to 2d array
pixels_2d = parse_to_2d(pixels, Nx, Ny, N_gst)
# define geometry at each pixel
geo = 0
for i, row in enumerate(pixels_2d):
for j, pixel in enumerate(row):
if pixel == 1:
geo += td.Box(
center=(l + Px / 2 + Px * j, (Ny - 1) / 2 * Py - Py * i, 0), size=(Px, Py, t)
)
# add input and output waveguides
geo += td.Box(center=(0, 0, 0), size=(2 * l, w, t))
geo += td.Box(center=(l + Nx * Px + l, 0, 0), size=(2 * l, w, t))
# define structure from geometry
design = td.Structure(geometry=geo, medium=Si)
return design
To test the functions, we define a parameter vector of all 1s. This should define a design that is all silicon besides the 3 by 3 pixels at the center for GST. The plot confirms that it's working correctly.
pixels = np.ones((Ny * Nx - N_gst * N_gst) // 2 + Nx // 2 - 1)
design = create_design(pixels)
design.plot(z=0)
plt.show()
Now we can define a function to create a Simulation given a parameter vector and the state of the GST.
def make_sim(pixels, gst_state):
# create silicon structures
design = create_design(pixels)
# add a mode source as excitation
mode_spec = td.ModeSpec(num_modes=1, target_neff=3.5)
mode_source = td.ModeSource(
center=(l / 2, 0, 0),
size=(0, 4 * w, 6 * t),
source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
direction="+",
mode_spec=mode_spec,
mode_index=0,
)
# add a mode monitor to measure transmission at the output waveguide
mode_monitor = td.ModeMonitor(
center=(Nx * Px + 1.5 * l, 0, 0),
size=mode_source.size,
freqs=freqs,
mode_spec=mode_spec,
name="mode",
)
# create the GST pixels
if gst_state == "c":
gst = td.Structure(
geometry=td.Box(center=(l + Px * Nx / 2, 0, 0), size=(N_gst * Px, N_gst * Py, t)),
medium=cGST,
)
elif gst_state == "a":
gst = td.Structure(
geometry=td.Box(center=(l + Px * Nx / 2, 0, 0), size=(N_gst * Px, N_gst * Py, t)),
medium=aGST,
)
else:
raise ValueError("gst_state can only be 'a' or 'c'.")
# define simulation domain size
Lx = 2 * l + Nx * Px
Ly = Ny * Py + 2 * buffer
Lz = 10 * t
# simulation domain box
sim_box = td.Box.from_bounds(rmin=(0, -Ly / 2, -Lz / 2), rmax=(Lx, Ly / 2, Lz / 2))
run_time = 1e-12 # simulation run time
# define simulation
sim = td.Simulation(
center=sim_box.center,
size=sim_box.size,
grid_spec=td.GridSpec.auto(min_steps_per_wvl=15, wavelength=lda0),
structures=[design, gst],
sources=[mode_source],
monitors=[mode_monitor],
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
medium=SiO2,
symmetry=(0, -1, 1),
)
return sim
Visualize the simulation in 3D to check if the setup is correct.
sim = make_sim(pixels, "c")
sim.plot_3d()
Optimization¶
The objective function is defined as the contrast between the transmissions when GST is in the amorphous and crystalline states. To evaluate it, we need to define two simulations with the same parameter vector but a different GST state, run the simulation in parallel, extract the transmission, and calculate the difference.
def objective_function(pixels):
# define two simulations
sim_a = make_sim(pixels, "a")
sim_c = make_sim(pixels, "c")
sims = {"a": sim_a, "c": sim_c}
# define and submit a batch
batch = web.Batch(simulations=sims, verbose=False)
batch_results = batch.run(path_dir="data")
# calculate transmission
sim_data_a = batch_results["a"]
sim_data_c = batch_results["c"]
t_a = sim_data_a["mode"].amps.sel(direction="+").values
T_a = np.mean(np.abs(t_a) ** 2)
t_c = sim_data_c["mode"].amps.sel(direction="+").values
T_c = np.mean(np.abs(t_c) ** 2)
return T_a - T_c
We start the optimization using the parameter vector of all 1s. The objective function is evaluated to be just above 0.8.
initial_obj = objective_function(pixels)
print(f"Initial objective function is {initial_obj:.3f}.")
Initial objective function is 0.812.
Since the DBS algorithm is pretty straightforward, we don't need to make sure of any external libraries. We can define the DBS optimizer directly. In the optimization, we flip the value of elements in the parameter vector one by one. If the flip leads to an improved objective function, we keep the flipped element. If no improvement is found, we revert back to the original value. One iteration ends when all elements of the parameter vector are flipped once and this can be done for multiple iterations. We can implement different termination conditions. For example, if no improvement is found or the improvement is smaller than a certain threshold between consecutive iterations, the optimization stops. Or if the objective function reaches a certain threshold value, the optimization stops. In this example, we simply monitor the improvement and terminate the optimization if no improvement is found within one iteration.
def direct_binary_search(pixels, iterations):
num_of_eval = 0 # number of objective function evaluations
best_score = [initial_obj] # best objective function values
print(f"The initial objective function is {initial_obj:.3f}.")
for i in range(iterations):
print(f"Iteration {i+1} starts.")
improvement = False # keep track of if there is improvement within one iteration
for idx in range(len(pixels)):
# flip the current element
pixels[idx] = 1 - pixels[idx]
# Evaluate the new solution
current_score = objective_function(pixels)
num_of_eval += 1
# if better, update best solution, else revert change
if current_score > best_score[-1]:
best_score.append(current_score)
improvement = True
print(
f"A best objective of {current_score:.3f} found after {num_of_eval} evaluation(s)."
)
else:
best_score.append(best_score[-1])
pixels[idx] = 1 - pixels[idx] # revert change
if not improvement:
break
return pixels, best_score
For simplicity, we only demonstrate the optimization for 1 iteration in this case. We tested more iterations and found that the second iteration only improved the performance very mildly. The third iteration showed no improvement at all. Within the first iteration, we improve the objective function by about 10%, which is quite significant.
pixel_opt, best_score = direct_binary_search(pixels, 1)
The initial objective function is 0.812. Iteration 1 starts. A best objective of 0.832 found after 1 evaluation(s). A best objective of 0.859 found after 2 evaluation(s). A best objective of 0.862 found after 6 evaluation(s). A best objective of 0.870 found after 7 evaluation(s). A best objective of 0.871 found after 8 evaluation(s). A best objective of 0.882 found after 9 evaluation(s). A best objective of 0.899 found after 10 evaluation(s). A best objective of 0.915 found after 11 evaluation(s). A best objective of 0.918 found after 13 evaluation(s). A best objective of 0.922 found after 18 evaluation(s). A best objective of 0.924 found after 24 evaluation(s). A best objective of 0.927 found after 26 evaluation(s).
Plot the best objective function as a function of the number of evaluations.
plt.plot(best_score, c="red", linewidth=2)
plt.xlabel("Number of evaluations")
plt.ylabel("Best objective function")
plt.show()
Final Design¶
After the optimization, we can inspect the final design more closely. We rerun the simulations with an added FieldMonitor so we can visualize the switching behavior better. We didn't have a FieldMonitor during the optimization to minimize unnecessary data download.
sim_a_opt = make_sim(pixels, "a")
sim_c_opt = make_sim(pixels, "c")
# add a field monitor to visualize field distribution at z=0
field_monitor = td.FieldMonitor(
center=(0, 0, 0), size=(td.inf, td.inf, 0), freqs=[freq0], name="field"
)
sim_a_opt = sim_a_opt.copy(update={"monitors": sim_a_opt.monitors + (field_monitor,)})
sim_c_opt = sim_c_opt.copy(update={"monitors": sim_c_opt.monitors + (field_monitor,)})
sims = {"a": sim_a_opt, "c": sim_c_opt}
batch_opt = web.Batch(simulations=sims, verbose=False)
batch_opt_results = batch_opt.run(path_dir="data")
Visualize the final design.
sim_a_opt.plot(z=0, monitor_alpha=0)
plt.show()
Plot the transmission spectra for the final design with GST at two states. A clear switching with >20 dB contrast is confirmed.
sim_data_a = batch_opt_results["a"]
sim_data_c = batch_opt_results["c"]
t_a = sim_data_a["mode"].amps.sel(direction="+").values
T_a = np.abs(t_a) ** 2
t_c = sim_data_c["mode"].amps.sel(direction="+").values
T_c = np.abs(t_c) ** 2
plt.plot(ldas, 10 * np.log10(T_a), c="blue", linewidth=2, label="On")
plt.plot(ldas, 10 * np.log10(T_c), c="red", linewidth=2, label="Off")
plt.xlabel("Wavelength (µm)")
plt.ylabel("Transmission (dB)")
plt.legend()
plt.show()
Finally, plot the electric field norms to see the energy flow. When the GST is in the amorphous state (on state), the waveguide mode is transmitted with little loss. When the GST is in the crystalline state (off state), the transmission is cut off nearly completely. The result confirms that the optimized design functions very well as an optical switch while maintaining an ultra-compact footprint.
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 2), tight_layout=True)
sim_data_a.plot_field(
field_monitor_name="field", field_name="E", val="abs", ax=ax1, eps_alpha=0.4, vmin=0, vmax=70
)
sim_data_c.plot_field(
field_monitor_name="field", field_name="E", val="abs", ax=ax2, eps_alpha=0.4, vmin=0, vmax=70
)
ax1.set_title("On state")
ax1.set_ylim(-1, 1)
ax2.set_title("Off state")
ax2.set_ylim(-1, 1)
plt.show()
