Note: the cost of running the entire notebook is about 10 FlexCredits.
Dielectric waveguides made of materials like silicon and other semiconductors have been well-known in integrated photonics to guide light with low loss. Plasmonic waveguides made of metals like silver and gold are another promising platform for nanoscale photonic applications relying on the use of surface plasmon polaritons (SPPs) to strongly confine and guide light. This tight field confinement leads to enhanced light-matter interactions, opening up applications in sensing, nanolasing, and quantum optics. Plasmonic waveguides such as metal-insulator-metal, insulator-metal-insulator, and metallic nanowires exhibit deep-subwavelength optical modes but also suffer from intrinsic losses due to metal absorption. Nonetheless, the combination of strong field confinement and sensitivity to surface conditions provides the potential for integrated nanophotonic devices and sensors based on plasmonic waveguides.
In this notebook, we model a plasmonic sensor for CO$_2$ concentration detection. The sensor is based on a metal-insulator-metal (MIM) waveguide coupled to a square ring cavity. The square ring cavity is filled with a polymer called polyhexamethylene biguanide (PHMB), which is sensitive to CO$_2$. When CO$_2$ is present, it causes a change in the refractive index of the PHMB, which shifts the resonant wavelength of the sensor. This shift can be used to detect and quantify the CO$_2$ concentration over a range of 0 - 524 ppm, which covers atmospheric concentrations relevant to climate change monitoring. The design is based on the work S. N. Khonina, N. L. Kazanskiy, M. A. Butt, A. Kaźmierczak, and R. Piramidowicz, "Plasmonic sensor based on metal-insulator-metal waveguide square ring cavity filled with functional material for the detection of CO2 gas," Opt. Express 29, 16584-16594 (2021) DOI:10.1364/OE.423141.
import numpy as np
import matplotlib.pyplot as plt
import tidy3d as td
import tidy3d.web as web
from tidy3d.plugins.mode import ModeSolver
Simulation Setup¶
Define the simulation wavelength range to be 1.1 μm to 1.6 μm.
lda0 = 1.35 # central wavelength
freq0 = td.C_0 / lda0 # central frequency
ldas = np.linspace(1.1, 1.6, 301) # wavelength range
freqs = td.C_0 / ldas # frequency range
fwidth = 0.5 * (np.max(freqs) - np.min(freqs)) # width of the source frequency
Four materials will be used in the simulation: gold, SiO$_2$, air, and the polymer PHMB. Gold is given by the Drude model with the parameters taken from the reference. SiO$_2$ is directly defined from the material library. Air has a constant permittivity of 1. PHMB has a refractive index sensitive to the CO$_2$ concentration, which will be investigated in a parameter sweep, so we will define the PHMB medium later.
# define gold medium
au = td.Drude(eps_inf=9.0685, coeffs=[(135.44e14 / (2 * np.pi), 1.15e14 / (2 * np.pi))])
# define sio2 medium
sio2 = td.material_library["SiO2"]["Palik_Lossless"]
# define air medium
air = td.Medium(permittivity=1)
Define the geometric parameters of the sensor. These values have been optimized in the reference. If further optimization is needed, parameter sweeps or adjoint optimization can be applied.
L = 0.6 # side length of the square cavity
g = 0.01 # gap size between the bus waveguide and the square cavity
W1 = 0.09 # width of the square cavity waveguide
W = 0.06 # width of the bus waveguide
H = 0.7 # thickness of the gold layer
inf_eff = 1e2 # effective infinity of the simulation
Define the structures of the sensor, which consists of the SiO$_2$ substrate, a gold layer, a bus waveguide filled with air, and a square ring cavity. The ring cavity is filled with PHMB. Therefore, we will define an outer square with PHMB and an inner square with gold. The outer square will be defined later with the PHMB medium.
In addition, we define a mesh override structure around the bus waveguide and square cavity such that the grid will be refined in this region. Since we want to accurately capture the resonance of the device, using a fine grid is important. The refractive index of the mesh override structure is set to 20.
# define the substrate structure
substrate = td.Structure(
geometry=td.Box.from_bounds(rmin=(-inf_eff, -inf_eff, -inf_eff), rmax=(inf_eff, inf_eff, 0)),
medium=sio2,
)
# define the gold layer structure on top of the substrate
au_layer = td.Structure(
geometry=td.Box.from_bounds(rmin=(-inf_eff, -inf_eff, 0), rmax=(inf_eff, inf_eff, H)), medium=au
)
# define the air-filled bus waveguide structure
bus_waveguide = td.Structure(
geometry=td.Box.from_bounds(rmin=(-W / 2, -inf_eff, 0), rmax=(W / 2, inf_eff, H)), medium=air
)
# define the inner square cavity structure
ring_cavity_inner = td.Structure(
geometry=td.Box(center=(W / 2 + g + L / 2, 0, H / 2), size=(L - 2 * W1, L - 2 * W1, H)),
medium=au,
)
# define the region where grid will be refined
refine_structure = td.Structure(
geometry=td.Box.from_bounds(rmin=(-W / 2, -L / 2, 0), rmax=(W / 2 + g + L, L / 2, H)),
medium=td.Medium(permittivity=20**2),
)
Define simulation domain size, run time, source, and monitors. Since the plasmonic waveguide sensor has a small footprint and the field is mostly confined inside the waveguide, the simulation domain only needs to be slightly larger than the device itself.
To excite the bus waveguide, a ModeSource is used, similar to simulating a dielectric waveguide. A ModeMonitor is defined at the other end of the bus waveguide to compute transmission. To visualize the field distribution, we also add a FieldMonitor at z=H/2. This monitor will record the fields at two frequencies, the minimum frequency and the central frequency. The minimum frequency is off-resonant and the central frequency is close to being on-resonant.
buffer = lda0 / 2 # buffer length in each direction
# define simulation domain size
Lx = L + buffer
Ly = L + g + W + buffer
Lz = H + buffer
run_time = 1.5e-13 # simulation run time
# add a mode source as excitation
mode_spec = td.ModeSpec(num_modes=1, target_neff=2)
mode_source = td.ModeSource(
center=(0, -L / 2 - lda0 / 4, H / 2),
size=(4 * W, 0, 3 * H),
source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
direction="+",
mode_spec=mode_spec,
mode_index=0,
num_freqs=5,
)
# add a mode monitor to measure transmission at the output waveguide
mode_monitor = td.ModeMonitor(
center=(0, L / 2 + lda0 / 4, H / 2),
size=mode_source.size,
freqs=freqs,
mode_spec=mode_spec,
name="mode",
)
# add a field monitor to visualize field distribution at z=H/2
field_monitor = td.FieldMonitor(
center=(0, 0, H / 2), size=(td.inf, td.inf, 0), freqs=[np.min(freqs), freq0], name="field"
)
[08:59:14] WARNING: Default value for the field monitor monitor.py:261 'colocate' setting has changed to 'True' in Tidy3D 2.4.0. All field components will be colocated to the grid boundaries. Set to 'False' to get the raw fields on the Yee grid instead.
To facilitate a parameter sweep to simulate different CO$_2$ concentrations, we will define a function make_sim(n_phmb), which returns a Tidy3D Simulation given the refractive index of PHMB. Since the simulation domain is small, if we use the PML boundary condition, we will receive warnings about structures being too close to the boundary. To avoid this warning, we use the Absorber boundary condition, which will work perfectly fine in this case. PML should be placed sufficiently far from any structures since evanescent field leaking into PML could cause the simulation to diverge. Absorber, on the other hand, does not have this concern.
def make_sim(n_phmb):
# define the phmb medium
phmb = td.Medium(permittivity=n_phmb**2)
# define the outer square cavity structure
ring_cavity_outer = td.Structure(
geometry=td.Box(center=(W / 2 + g + L / 2, 0, H / 2), size=(L, L, H)), medium=phmb
)
# define simulation
sim = td.Simulation(
center=(W / 2 + g + L / 2, 0, H / 2),
size=(Lx, Ly, Lz),
grid_spec=td.GridSpec.auto(
min_steps_per_wvl=30, wavelength=lda0, override_structures=[refine_structure]
),
structures=[substrate, au_layer, bus_waveguide, ring_cavity_outer, ring_cavity_inner],
sources=[mode_source],
monitors=[mode_monitor, field_monitor],
run_time=run_time,
boundary_spec=td.BoundarySpec.all_sides(boundary=td.Absorber()),
)
return sim
To test our make_sim function, we define a simulation and visualize it.
sim = make_sim(1.55)
sim.plot(z=H / 2)
plt.show()
We further visualize the grid around the cavity region to ensure it's sufficiently fine.
ax = sim.plot(z=H / 2)
sim.plot_grid(z=H / 2, ax=ax)
ax.set_xlim(-W / 2, W / 2 + g + W1 + L)
ax.set_ylim(-L / 2, L / 2)
ax.set_aspect("auto")
plt.show()
Before running the simulations, we can use the ModeSolver plugin to visualize the source mode to ensure we are launching the plasmonic waveguide mode.
mode_solver = ModeSolver(
simulation=sim,
plane=td.Box(center=mode_source.center, size=mode_source.size),
mode_spec=mode_spec,
freqs=[freq0],
)
mode_data = mode_solver.solve()
[08:59:16] WARNING: Use the remote mode solver with subpixel mode_solver.py:154 averaging for better accuracy through 'tidy3d.plugins.mode.web.run(...)'.
After running the solver, we can visualize the mode field. Indeed we see the plasmonic waveguide mode with the field concentrated inside the bus waveguide.
mode_data.intensity.plot(x="x", y="z", cmap="magma", vmin=0, vmax=3e4)
plt.show()
Lastly, we can estimate the cost of running this single simulation, which will give us a better idea of the cost of running a batch of simulations.
job = web.Job(simulation=sim, task_name="plasmonic_sensor")
estimated_cost = web.estimate_cost(job.task_id)
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- dc128cb3-c2f8-4bf2-9440-9816e5adeb0av1'.
Output()
The estimated maximum cost is 2.633 Flex Credits.
Running the Parameter Sweep¶
Now that everything is set up correctly, we are ready to run a parameter sweep to investigate the transmission spectra at different PHMB refractive index values corresponding to different CO$_2$ concentrations. For the concentration in the range of 0 ppm to 524 ppm, the refractive index varies linearly from 1.55 to 1.48.
# refractive index values to be simulated
n_phmb_array = np.array([1.48, 1.51, 1.53, 1.55])
# define a simulation batch
sims = {f"n_phmb={n_phmb:.2f}": make_sim(n_phmb) for n_phmb in n_phmb_array}
batch = web.Batch(simulations=sims)
# run the batch
batch_results = batch.run(path_dir="data")
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 586fb1e0-6257-47fc-9647-1834a9d21b09v1'.
Output()
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 2b6f5d9f-9c8d-4471-bffa-e727ef6d9d8fv1'.
Output()
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- 67c1ce7b-8a42-400b-92e8-43b03731d64ev1'.
Output()
View task using web UI at webapi.py:190 'https://tidy3d.simulation.cloud/workbench?taskId=fdve- e7236603-b480-4138-9a29-0a54a5d9d9dav1'.
Output()
[08:59:46] Started working on Batch. container.py:475
[08:59:49] Maximum FlexCredit cost: 10.531 for the whole batch. container.py:479 Use 'Batch.real_cost()' to get the billed FlexCredit cost after the Batch has completed.
Output()
[09:02:35] Batch complete. container.py:522
After running the parameter sweep, we can check the actual FlexCredit cost.
real_cost = batch.real_cost()
The actual FlexCredit cost is 9.899322401520688.
Result Visualization¶
After the simulations are complete, we can plot the transmission spectra at different PHMB refractive index values. We see a significant blueshift of the resonant wavelength as the refractive index decreases from 1.55 to 1.48. Therefore, we can use the resonance wavelength to determine the concentration of CO$_2$.
for n_phmb in n_phmb_array:
sim_data = batch_results[f"n_phmb={n_phmb:.2f}"]
amp = sim_data["mode"].amps.sel(mode_index=0, direction="+")
T = np.abs(amp) ** 2
plt.plot(ldas, 10 * np.log10(T), label=f"n={n_phmb}")
plt.legend()
plt.xlim(np.min(ldas), np.max(ldas))
plt.ylim(-35, -3)
plt.xlabel("Wavelength ($\mu m$)")
plt.ylabel("Transmission (dB)")
plt.show()
Output()
Output()
Output()
Output()
Finally, we plot the field ($H_z$) distributions on- and off-resonance. When it's off-resonance, we see a high transmission of the plasmonic waveguide mode while a low transmission when it's on resonance, consistent with the spectra above.
sim_data = batch_results[f"n_phmb=1.51"]
fig, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(8, 4))
sim_data.plot_field(
field_monitor_name="field",
field_name="Hz",
val="real",
ax=ax1,
f=np.min(freqs),
vmin=-0.3,
vmax=0.3,
)
sim_data.plot_field(
field_monitor_name="field", field_name="Hz", val="real", ax=ax2, f=freq0, vmin=-0.3, vmax=0.3
)
plt.show()
To fully characterize the performance of a refractive index sensor, one can perform further quantitative analysis on the resonance wavelength shift and compute quantities like the sensitivity. For the sake of simplicity, we will not delve into it in this notebook.
