TIDY3D
LEARNING CENTER

Lectures
Tutorials
Applications

Example Library / Mode Solver / Mode confinement factor of slot and strip waveguides

Mode confinement factor of slot and strip waveguides

Sign up now to simulate .ipynb

This tutorial computes the mode confinement factor of dielectric waveguides with the Tidy3D ModeSimulation, for silicon-on-insulator (SOI) strip and slot geometries. The confinement factor measures how much of the guided mode overlaps a chosen region, here the low-index air cladding, which is a key quantity for sensing and modulator design. The calculation follows Kita et al., Optica 5, 1046 (2018), using equation 3 and the geometries of figure 1.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import tidy3d as td
import tidy3d.web as web

Confinement Factor

Following Kita et al., we compute two confinement factors. The external (bulk) factor is the fraction of modal energy in the low-index air cladding (equation 3),

\[\Gamma_\text{bulk} = \frac{n_g}{n_{\text{air}}}\,\frac{\int_{\text{air}}\epsilon|E|^2\,dA}{\int\epsilon|E|^2\,dA},\]

while the surface factor uses the same expression but integrates only over an 8 nm thick air shell along the air/solid interface, capturing the field available to surface-bound analytes. Here \(n_g\) is the group index. The helper below returns both, given masks for the silicon core and for the surface shell.

SHELL = 0.008   # 8 nm air shell along the air/solid interface (Kita surface region)


def box_dist(X, Y, x1, x2, y1, y2):
    dx = np.maximum(np.maximum(x1 - X, X - x2), 0)
    dy = np.maximum(np.maximum(y1 - Y, Y - y2), 0)
    return np.hypot(dx, dy)


def confinement(m, i, is_core, is_surf):
    f = lambda c: getattr(m, c).sel(mode_index=i, f=freq, method="nearest").squeeze().transpose("x", "y").values
    E2 = np.abs(f("Ex"))**2 + np.abs(f("Ey"))**2 + np.abs(f("Ez"))**2
    x, y = m.Ex.x.values, m.Ex.y.values
    X, Y = np.meshgrid(x, y, indexing="ij")
    I = lambda mask: np.trapezoid(np.trapezoid(np.where(mask, E2, 0), y, 1), x, 0)
    core = I(is_core(X, Y))
    box = I(Y <= 0)                        # SiO2 substrate
    air = I(True) - core - box
    ng = float(m.n_group.sel(mode_index=i, f=freq, method="nearest").real)
    W = 3.48**2 * core + 1.44**2 * box + air   # n_air = 1
    return ng * air / W, ng * I(is_surf(X, Y)) / W   # bulk, surface

Setup

We define the material and geometry parameters shared by both waveguide families: a 220 nm silicon core (\(n = 3.48\)) on a silica (\(\text{SiO}_2\), \(n = 1.44\)) substrate at a wavelength of 1.55 µm. Air is the background Medium (\(n = 1\)), so the slot gap needs no explicit structure.

h = 0.22                      # Si thickness (um)
wl = 1.55                     # wavelength (um)
freq = td.C_0 / wl
si = td.Medium(permittivity=3.48**2)
sio2 = td.Medium(permittivity=1.44**2)

Mode Selection

The fundamental TE and TM modes are picked with a ModeSortSpec that orders the computed modes by the field fraction inside the waveguide box (filter_order="over" keeps the TE branch first, "under" the TM branch).

def spec(box, order):
    return td.ModeSortSpec(filter_key="TE_fraction", filter_reference=0.5, filter_order=order,
                           sort_key="fill_fraction_box", sort_order="descending", bounding_box=box)


def sim(structs, wg_box, mesh_box):
    return td.ModeSimulation(
        size=(2, 3, 1),
        structures=structs,
        grid_spec=td.GridSpec.auto(wavelength=wl, min_steps_per_wvl=50,
            override_structures=[td.MeshOverrideStructure(geometry=mesh_box, dl=(1e-3, 1e-3, 1e-2))]),
        boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),   # absorb substrate/handle radiation
        mode_spec=td.ModeSpec(num_modes=10, target_neff=2.0, num_pml=(9, 9), precision="double",
                              group_index_step=True, sort_spec=spec(wg_box, "over")),   # rank by field fraction inside the Si waveguide
        freqs=[freq],
        plane=td.Box(center=(0, h / 2, 0), size=(2, 3, 0)),
    )


def extract_confinement(results, specs, xname):
    rows = []
    for k, sd in results.items():
        v = specs[k]
        te = sd.sort_modes(spec(v["wg"], "over"))                    # most waveguide-confined TE
        tm = sd.sort_modes(spec(v["wg"], "under"))                   # most waveguide-confined TM
        te_b, te_s = confinement(te.modes, 0, v["core"], v["surf"])
        tm_b, tm_s = confinement(tm.modes, 0, v["core"], v["surf"])
        rows.append({xname: v["x"], "width": v["width"],
                     "TE_bulk": te_b, "TE_surf": te_s, "TM_bulk": tm_b, "TM_surf": tm_s})
    return pd.DataFrame(rows).sort_values([xname, "width"])

Slot Waveguide

Two silicon rails separated by an air slot on a silica substrate. We sweep the slot width from 40 nm to 160 nm for three total waveguide widths (450, 550, and 650 nm). Before running, we preview one ModeSimulation per total width.

# --- Slot: two Si rails on a SiO2 substrate; air gap = background ---
widths = {"450": 0.45, "550": 0.55, "650": 0.65}
slot_specs = {}
for wl_, w in widths.items():
    for s in np.round(np.arange(0.04, 0.161, 0.01), 3):
        s, rail = float(s), (w - float(s)) / 2
        wg = td.Box(center=(0, h / 2, 0), size=(w, h, 0))       # whole waveguide, for mode selection
        mesh = td.Box(center=(0, h / 2, 0), size=(w + 4 * SHELL, h + 4 * SHELL, 0))   # guide + 8 nm surface shell
        structs = [
            td.Structure(geometry=td.Box(center=(0, -5, 0), size=(td.inf, 10, td.inf)), medium=sio2),   # semi-infinite substrate
            td.Structure(geometry=td.Box(center=(-(s / 2 + rail / 2), h / 2, 0), size=(rail, h, td.inf)), medium=si),
            td.Structure(geometry=td.Box(center=(s / 2 + rail / 2, h / 2, 0), size=(rail, h, td.inf)), medium=si),
        ]
        core = lambda X, Y, s=s, rail=rail: (Y > 0) & (Y <= h) & (
            ((X >= -(s / 2 + rail)) & (X <= -s / 2)) | ((X >= s / 2) & (X <= s / 2 + rail)))
        surf = lambda X, Y, s=s, rail=rail, core=core: (~(core(X, Y, s, rail) | (Y <= 0))) & (np.minimum(np.minimum(
            box_dist(X, Y, -(s / 2 + rail), -s / 2, 0, h),
            box_dist(X, Y, s / 2, s / 2 + rail, 0, h)), np.maximum(Y, 0.0)) <= SHELL)
        slot_specs[f"w{wl_}_s{s:.3f}"] = dict(sim=sim(structs, wg, mesh), wg=wg, core=core, surf=surf, x=s, width=w)

# preview one ModeSimulation per total width before running the sweep
fig, axes = plt.subplots(1, len(widths), figsize=(12, 3.5), tight_layout=True)
for ax, (wl_, w) in zip(axes, widths.items()):
    slot_specs[f"w{wl_}_s0.080"]["sim"].plot(z=0, ax=ax)
    ax.set_title(f"total width = {wl_} nm")
plt.show()

slot_results = web.Batch(simulations={k: v["sim"] for k, v in slot_specs.items()}).run(path_dir="data_slot")
slot = extract_confinement(slot_results, slot_specs, "slot")







07:02:55 UTC Started working on Batch containing 39 tasks.                      







07:03:26 UTC Maximum FlexCredit cost: 3.380 for the whole batch.                







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












07:06:05 UTC Batch complete.                                                    

















Strip Waveguide


A single silicon core on a silica substrate, with the width swept from 300 nm to 500 nm. We preview one ModeSimulation per width before running the sweep.






# --- Strip: one Si core on a SiO2 substrate ---
strip_specs = {}
for w in np.round(np.linspace(0.30, 0.50, 5), 3):
    w = float(w)
    wg = td.Box(center=(0, h / 2, 0), size=(w, h, 0))       # core, for mode selection
    mesh = td.Box(center=(0, h / 2, 0), size=(w + 4 * SHELL, h + 4 * SHELL, 0))   # core + 8 nm surface shell
    structs = [
        td.Structure(geometry=td.Box(center=(0, -5, 0), size=(td.inf, 10, td.inf)), medium=sio2),   # semi-infinite substrate
        td.Structure(geometry=td.Box(center=(0, h / 2, 0), size=(w, h, td.inf)), medium=si),
    ]
    core = lambda X, Y, w=w: (Y > 0) & (Y <= h) & (np.abs(X) <= w / 2)
    surf = lambda X, Y, w=w, core=core: (~(core(X, Y, w) | (Y <= 0))) & (
        np.minimum(box_dist(X, Y, -w / 2, w / 2, 0, h), np.maximum(Y, 0.0)) <= SHELL)
    strip_specs[f"w{w:.3f}"] = dict(sim=sim(structs, wg, mesh), wg=wg, core=core, surf=surf, x=w, width=w)

# preview one ModeSimulation per width before running the sweep
fig, axes = plt.subplots(1, len(strip_specs), figsize=(15, 3), tight_layout=True)
for ax, (k, v) in zip(axes, strip_specs.items()):
    v["sim"].plot(z=0, ax=ax)
    ax.set_title(f"width = {int(round(v['width'] * 1000))} nm")
plt.show()

strip_results = web.Batch(simulations={k: v["sim"] for k, v in strip_specs.items()}).run(path_dir="data_strip")
strip = extract_confinement(strip_results, strip_specs, "width")
































07:13:54 UTC Started working on Batch containing 5 tasks.                       







07:13:58 UTC Maximum FlexCredit cost: 0.343 for the whole batch.                







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












07:16:41 UTC Batch complete.                                                    

















Results


We plot the bulk and surface confinement factors: the strip against its width, and the slot against the slot width, with a different marker for each total waveguide width.






# --- Bulk and surface confinement factors ---
mk = {"450": "^", "550": "s", "650": "x"}
fig, ax = plt.subplots(2, 2, figsize=(13, 9))

# Strip (left column): TE/TM, bulk on top and surface below
for br, c in [("TE", "tab:red"), ("TM", "tab:blue")]:
    ax[0, 0].plot(strip["width"] * 1e3, strip[f"{br}_bulk"], "o-", color=c, ms=7, label=br)
    ax[1, 0].plot(strip["width"] * 1e3, strip[f"{br}_surf"], "o-", color=c, ms=7)
ax[0, 0].set(title="Strip", ylabel="Bulk Γ", ylim=(0, 1))
ax[1, 0].set(xlabel="Strip width (nm)", ylabel="Surface Γ", ylim=(0, 0.08))
ax[0, 0].legend()

# Slot (right column): one marker per total width, color per polarization
for wl_, w in widths.items():
    d = slot[slot.width == w]
    for br, c in [("TE", "tab:red"), ("TM", "tab:blue")]:
        lbl = f"width = {wl_} nm" if br == "TE" else None
        ax[0, 1].plot(d["slot"] * 1e3, d[f"{br}_bulk"], marker=mk[wl_], color=c, lw=1.2, ms=7, label=lbl)
        ax[1, 1].plot(d["slot"] * 1e3, d[f"{br}_surf"], marker=mk[wl_], color=c, lw=1.2, ms=7)
ax[0, 1].set(title="Slot", ylabel="Bulk Γ", ylim=(0, 1))
ax[1, 1].set(xlabel="Slot width (nm)", ylabel="Surface Γ", ylim=(0, 0.3))
ax[0, 1].legend(fontsize=8)

for a in ax.flat:
    a.grid(ls="--", alpha=0.4)
fig.tight_layout()
plt.show()





















Both the bulk and surface confinement factors reproduce the trends of figure 1 in Kita et al. (2018): the slot maximizes the TE overlap with the low-index gap, while for the strip the TM mode extends further into the cladding than the TE mode.






Mode Field Profiles


Electric field magnitude \(|E|\) of the fundamental TE and TM modes at the first and last point of each sweep (strip 300 and 500 nm; slot 450 nm/40 nm and 650 nm/160 nm).






# --- |E| of the fundamental TE and TM modes at the first and last sweep points ---
cases = {"strip 300 nm": strip_specs["w0.300"], "strip 500 nm": strip_specs["w0.500"],
         "slot 450/40 nm": slot_specs["w450_s0.040"], "slot 650/160 nm": slot_specs["w650_s0.160"]}
fdata = web.Batch(simulations={k: v["sim"] for k, v in cases.items()}).run(path_dir="data_fields")

fig, ax = plt.subplots(len(cases), 2, figsize=(10, 3.4 * len(cases)), tight_layout=True)
for i, (k, v) in enumerate(cases.items()):
    sd = fdata[k]
    for j, order in enumerate(["over", "under"]):
        m = (sd if order == "over" else sd.sort_modes(spec(v["wg"], "under"))).modes
        E = (abs(m.Ex) ** 2 + abs(m.Ey) ** 2 + abs(m.Ez) ** 2) ** 0.5
        E.sel(mode_index=0).isel(f=0).squeeze().T.plot(ax=ax[i, j], cmap="magma", add_colorbar=False)
        ax[i, j].set_title(f"{k}, {'TE' if j == 0 else 'TM'}")
        ax[i, j].set_xlim(-1, 1)
        ax[i, j].set_ylim(-1.2, 1.2)
plt.show()


















07:17:01 UTC Started working on Batch containing 4 tasks.                       







07:17:04 UTC Maximum FlexCredit cost: 0.311 for the whole batch.                







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












07:18:53 UTC Batch complete.                                                    






























          

          

            Example library
          

          

          
          

          

            
            
              
              Previous Tutorial
            
            
            
            
              Next Tutorial