TIDY3D
LEARNING CENTER

Lectures
Tutorials
Applications

Example Library / Inverse Design / Compact TE modal demultiplexer: topology optimization with contour constraints

Compact TE modal demultiplexer: topology optimization with contour constraints

Sign up now to simulate .ipynb

This notebook reproduces the compact modal demultiplexer from §III-B of Wong et al., “Simplifying Photonic Topologies in Inverse Design With Contour Constraints” [1]. The device is a 5 × 4 µm² O-band TE modal demultiplexer on a 220 nm silicon-on-insulator (SOI) platform with 700 nm access waveguides.

The contour constraint is a foundry-agnostic topology-optimization (TO) technique that automatically identifies and erodes connected regions carrying low integrated field power — the “extraneous” features that inflate device complexity without contributing to performance. It removes the manual post-processing step of deleting non-functional features by hand, yielding simpler, more fabrication-friendly geometries with no penalty to insertion loss.

diagram

Routing function

Input mode Output port Output mode
TE0 Bottom TE0
TE1 Top TE0 (mode-converted from TE1)

Two designs compared over the O-band (1260–1360 nm):

  1. Standard TO demux — density-based TO with foundry constraints.
  2. Contour-constrained TO demux — the same, plus the contour constraint that gradually erodes connected regions with low integrated field power.

Optimization workflow (following the paper’s Fig. 2 milestones):

  • Phase 1 (staged TO) — unconstrained TO (iters 0–50) → ramped fabrication penalty (50–120) → full fabrication constraints (120–135). Shared between both branches.
  • Phase 2 (parallel branches) — from iteration 135, standard TO continued versus contour-constrained TO.
  • Phase 3 (seeded TO) — a fabrication-only polishing pass applied to both branches.

This example uses the Tidy3D autograd plugin for automatic-differentiation–based inverse design. If you are new to Tidy3D, we recommend starting with the tutorials and browsing the example library for more inverse design walkthroughs.

import autograd
import autograd.numpy as np
import matplotlib.pyplot as plt
import scipy.ndimage

import tidy3d as td
import tidy3d.web as web
from scipy.interpolate import RegularGridInterpolator
from tidy3d.plugins.autograd import (
    make_filter_and_project,
    make_erosion_dilation_penalty,
)
from tidy3d.plugins.invdes import AdamOptimizer

1. O-Band SOI Platform Parameters

Matches the paper: 220 nm SOI, 700 nm waveguide width, O-band (1260–1360 nm) on the GF Fotonix™ platform.

# O-band wavelength range
lda0 = 1.31  # central wavelength (μm)
freq0 = td.C_0 / lda0

lda_min = 1.26
lda_max = 1.36
freq_min = td.C_0 / lda_max
freq_max = td.C_0 / lda_min
fwidth = 0.5 * (freq_max - freq_min)

# Materials
n_si = 3.476
n_sio2 = 1.444
eps_si = n_si**2
si = td.Medium(permittivity=eps_si)
sio2 = td.Medium(permittivity=n_sio2**2)

# Geometry (paper values)
wg_width = 0.7  # 700 nm waveguides
thick = 0.22  # 220 nm SOI silicon thickness

# Frequency arrays
num_analysis_freqs = 81
ldas_analysis = np.linspace(lda_min, lda_max, num_analysis_freqs)
freqs_analysis = td.C_0 / ldas_analysis

num_opt_freqs = 5
ldas_opt = np.linspace(lda_min, lda_max, num_opt_freqs)
freqs_opt = list(td.C_0 / ldas_opt)

print("O-band SOI demux setup")
print(f"  Waveguide width: {wg_width*1000:.0f} nm")
print(f"  Si thickness:    {thick*1000:.0f} nm")
print(f"  Wavelength band: {lda_min:.2f}–{lda_max:.2f} μm")
print(f"  Optimization frequencies: {num_opt_freqs} pts")
O-band SOI demux setup
  Waveguide width: 700 nm
  Si thickness:    220 nm
  Wavelength band: 1.26–1.36 μm
  Optimization frequencies: 5 pts

2. Demultiplexer Geometry

  • 5 µm (x) × 4 µm (y) design region
  • 1 input waveguide on the left (centered at y=0)
  • 2 output waveguides on the right (top: y=+1.0 µm, bottom: y=−1.0 µm)
  • 1 µm of straight waveguide on each side beyond the design region

The 1.0 µm offset between output ports keeps them well-separated relative to the mode field diameter (≈0.7 µm) so we don’t get spurious evanescent coupling outside the design region.

# Design region — rectangular
design_size_x = 5.0  # μm
design_size_y = 4.0  # μm
wg_length = 1.0  # straight waveguide length each side
pixel_size = 0.025  # 25 nm
radius = 0.1  # 100 nm minimum feature size (filter radius)
min_steps_per_wvl = lda0 / (pixel_size * np.sqrt(eps_si))

num_pixels_x = int(design_size_x / pixel_size)
num_pixels_y = int(design_size_y / pixel_size)

# Output port y-offset
y_port = 1.0

print(f"Design region: {design_size_x} × {design_size_y} μm²")
print(
    f"Pixel grid:    {num_pixels_x} × {num_pixels_y}  ({num_pixels_x*num_pixels_y} voxels)"
)
print(f"Output port offset: ±{y_port} μm")

# Static waveguide structures
waveguide_in = td.Structure(
    geometry=td.Box(
        center=(-wg_length - design_size_x / 2, 0, 0),
        size=(2 * wg_length, wg_width, thick),
    ),
    medium=si,
)

waveguide_top = td.Structure(
    geometry=td.Box(
        center=(+wg_length + design_size_x / 2, +y_port, 0),
        size=(2 * wg_length, wg_width, thick),
    ),
    medium=si,
)

waveguide_bot = td.Structure(
    geometry=td.Box(
        center=(+wg_length + design_size_x / 2, -y_port, 0),
        size=(2 * wg_length, wg_width, thick),
    ),
    medium=si,
)

design_region_geometry = td.Box(
    center=(0, 0, 0),
    size=(design_size_x, design_size_y, thick),
)
Design region: 5.0 × 4.0 μm²
Pixel grid:    200 × 160  (32000 voxels)
Output port offset: ±1.0 μm

3. Sources and Monitors

To define the dual-mode FOM we run two simulations per gradient step:

  • Sim A: TE0 source on the input → measure TE0 amplitude at the bottom monitor
  • Sim B: TE1 source on the input → measure TE0 amplitude at the top monitor

Both monitors use num_modes=2 so we can examine cross-talk into the other output port if needed.

src_size_y = 4 * wg_width  # generous to fully capture both TE0 and TE1
src_size_z = 6 * thick

# Source x-position: just inside the input straight section
x_src = -design_size_x / 2 - wg_length + lda0 / 3
# Monitor x-position: just inside each output straight section
x_mnt = +design_size_x / 2 + wg_length - lda0 / 3

mode_spec_input = td.ModeSpec(num_modes=2, target_neff=n_si)
mode_spec_output = td.ModeSpec(num_modes=2, target_neff=n_si)


def make_source(mode_index: int) -> td.ModeSource:
    return td.ModeSource(
        center=(x_src, 0, 0),
        size=(0, src_size_y, src_size_z),
        source_time=td.GaussianPulse(freq0=freq0, fwidth=fwidth),
        direction="+",
        mode_spec=mode_spec_input,
        mode_index=mode_index,
    )


mode_source_te0 = make_source(0)
mode_source_te1 = make_source(1)


# Output monitors (one per port) — record both TE0 and TE1 at each
def make_monitor(y: float, name: str, freqs=None) -> td.ModeMonitor:
    return td.ModeMonitor(
        center=(x_mnt, y, 0),
        size=(0, src_size_y, src_size_z),
        freqs=list(freqs_opt) if freqs is None else list(freqs),
        mode_spec=mode_spec_output,
        name=name,
    )


mode_monitor_top = make_monitor(+y_port, "mode_top")
mode_monitor_bot = make_monitor(-y_port, "mode_bot")

# Optional 2D field monitor for visualization
field_monitor = td.FieldMonitor(
    center=(0, 0, 0),
    size=(td.inf, td.inf, 0),
    freqs=[freq0],
    name="field",
)

4. Topology Optimization Setup

Density-based TO with a conic filter (minimum feature size = 100 nm) and a hyperbolic-tangent projection. The design region holds the demultiplexer’s nontrivial topology.

filter_project_fn = make_filter_and_project(radius=radius, dl=pixel_size)


def get_density(params: np.ndarray, beta: float) -> np.ndarray:
    return filter_project_fn(params, beta=beta)


def get_design_region(params: np.ndarray, beta: float) -> td.Structure:
    density = get_density(params, beta=beta)
    eps_data = n_sio2**2 + (eps_si - n_sio2**2) * density
    return td.Structure.from_permittivity_array(
        eps_data=eps_data, geometry=design_region_geometry
    )


# Base simulation — the actual sources/monitors are swapped in per call
sim_base = td.Simulation(
    size=(
        2 * wg_length + design_size_x,
        2 * wg_length + design_size_y,
        thick + 2 * lda0,
    ),
    run_time=120 / fwidth,
    structures=[waveguide_in, waveguide_top, waveguide_bot],
    boundary_spec=td.BoundarySpec.all_sides(boundary=td.PML()),
    grid_spec=td.GridSpec.auto(
        min_steps_per_wvl=min_steps_per_wvl,
        override_structures=[
            td.MeshOverrideStructure(
                geometry=design_region_geometry, dl=3 * [pixel_size], enforce=True
            )
        ],
    ),
    medium=sio2,
    sources=[mode_source_te0],
    monitors=[mode_monitor_top, mode_monitor_bot],
)


def get_sim(
    params: np.ndarray,
    beta: float,
    input_mode: int,
    broadband: bool = False,
    with_field: bool = False,
) -> td.Simulation:
    """Build a simulation with the requested input mode and monitor frequencies."""
    design_region = get_design_region(params, beta=beta)
    src = mode_source_te0 if input_mode == 0 else mode_source_te1
    if broadband:
        m_top = make_monitor(+y_port, "mode_top", freqs=freqs_analysis)
        m_bot = make_monitor(-y_port, "mode_bot", freqs=freqs_analysis)
    else:
        m_top, m_bot = mode_monitor_top, mode_monitor_bot
    monitors = [m_top, m_bot]
    if with_field:
        monitors = monitors + [field_monitor]
    return sim_base.updated_copy(
        structures=list(sim_base.structures) + [design_region],
        sources=[src],
        monitors=monitors,
    )


# Initial parameters — uniform 0.5 (paper's standard starting point)
params0 = 0.5 * np.ones((num_pixels_x, num_pixels_y, 1))

print(f"Total design parameters: {params0.size}")
Total design parameters: 32000
# Visualize the initial geometry
sim_te0_init = get_sim(params0, beta=8.0, input_mode=0)

fig, ax = plt.subplots(figsize=(10, 5))
sim_te0_init.plot_eps(z=0.01, ax=ax)
ax.set_title("Initial demux geometry (top view)")
plt.tight_layout()
plt.show()

5. Dual-Mode Objective

For the modal demultiplexer the figure of merit is the average of two broadband transmission terms, one per input mode:

\[ \mathrm{FOM} = \tfrac{1}{2}\Big( \overline{|s^{\,(\text{TE}_0 \text{ in})}_{\text{bot},\,\text{TE}_0}|^2} \;+\; \overline{|s^{\,(\text{TE}_1 \text{ in})}_{\text{top},\,\text{TE}_0}|^2} \Big) \]

Each step requires two FDTD runs (one per input mode). Autograd flows through both and gradients are summed automatically.

_last_transmission = [None]  # (T_te0, T_te1) for live IL printing
_last_T_te0 = [None]
_last_T_te1 = [None]


def get_transmission(params: np.ndarray, beta: float) -> float:
    """Dual-simulation FOM: averages TE0→bot and TE1→top transmissions."""
    sim_te0 = get_sim(params, beta=beta, input_mode=0)
    sim_te1 = get_sim(params, beta=beta, input_mode=1)

    # Batch run for parallel cloud execution
    batch_data = web.run_async(
        {"te0_in": sim_te0, "te1_in": sim_te1},
        folder_name="contour constraint demux",
        verbose=False
    )

    amp_te0_to_bot = (
        batch_data["te0_in"]["mode_bot"].amps.sel(mode_index=0, direction="+").values
    )
    amp_te1_to_top = (
        batch_data["te1_in"]["mode_top"].amps.sel(mode_index=0, direction="+").values
    )

    T_te0 = np.mean(np.abs(amp_te0_to_bot) ** 2)
    T_te1 = np.mean(np.abs(amp_te1_to_top) ** 2)

    raw_te0 = T_te0._value if hasattr(T_te0, "_value") else T_te0
    raw_te1 = T_te1._value if hasattr(T_te1, "_value") else T_te1
    _last_T_te0[0] = float(raw_te0)
    _last_T_te1[0] = float(raw_te1)

    return 0.5 * (T_te0 + T_te1)


# Erosion-dilation fabrication penalty
over_under_etch = 0.1
penalty_fn = make_erosion_dilation_penalty(
    radius=over_under_etch, dl=pixel_size, beta=10.0
)


def get_penalty(params: np.ndarray, beta: float) -> float:
    density = get_density(params, beta=beta)
    return penalty_fn(density)


# Staged fab weight — set per Phase-1 iteration
_fab_weight = [1.0]


def objective(params: np.ndarray, beta: float) -> float:
    """Transmission − weighted fabrication penalty."""
    transmission = get_transmission(params, beta=beta)
    penalty = get_penalty(params, beta=beta)
    return transmission - _fab_weight[0] * penalty

6. Phase 1: Staged Standard Topology Optimization

Phase 1 follows the paper’s Fig. 2 schedule:

Iters Fab weight β Stage
0–49 0.0 8 → ~16 Unconstrained — fields not yet confined, extraneous contours form
50–119 0 → 1 (ramp) ~16 → ~28 Foundry constraints fade in
120–134 1.0 ~28 → 32 Full DRC

This is critical: without an unconstrained early stage the optimizer converges directly to a clean topology and there are no extraneous contours for the constraint to act on later.

num_steps_phase1 = 75
learning_rate = 0.02

beta_min = 8
beta_max = 32


def get_beta(step_num: int, num_steps: int) -> float:
    return beta_min + (beta_max - beta_min) * step_num / max(num_steps - 1, 1)


unconstrained_end = 25
ramp_end = 65


def get_fab_weight(step_num: int) -> float:
    if step_num < unconstrained_end:
        return 0.0
    if step_num < ramp_end:
        return (step_num - unconstrained_end) / (ramp_end - unconstrained_end)
    return 1.0


params = params0.copy()
optimizer_adam = AdamOptimizer.model_construct(learning_rate=learning_rate)
opt_state = optimizer_adam.initial_state(params)

val_grad_fn = autograd.value_and_grad(objective)

objective_history_p1 = []
il_history_p1 = []  # average IL across both inputs
il_te0_history_p1 = []
il_te1_history_p1 = []
param_history_p1 = [params.copy()]

print("=" * 70)
print("PHASE 1: Staged TO (unconstrained → ramp → full fab)")
print("=" * 70)

for i in range(num_steps_phase1):
    beta = get_beta(i, num_steps_phase1)
    _fab_weight[0] = get_fab_weight(i)
    stage = (
        "unconstrained"
        if i < unconstrained_end
        else ("ramping" if i < ramp_end else "full-fab")
    )

    value, gradient = val_grad_fn(params, beta=beta)

    params, opt_state = optimizer_adam.update(
        parameters=params, gradient=-gradient, state=opt_state
    )
    params = np.clip(params, 0, 1)

    il_te0 = -10 * np.log10(max(_last_T_te0[0], 1e-12))
    il_te1 = -10 * np.log10(max(_last_T_te1[0], 1e-12))
    il_avg = 0.5 * (il_te0 + il_te1)

    objective_history_p1.append(float(value))
    il_history_p1.append(il_avg)
    il_te0_history_p1.append(il_te0)
    il_te1_history_p1.append(il_te1)
    param_history_p1.append(params.copy())

    print(
        f"  Step {i+1:3d}/{num_steps_phase1}, β={beta:.1f}, fab_w={_fab_weight[0]:.2f} [{stage}]"
        f"  obj={value:.6f}  IL(TE0)={il_te0:.4f} dB  IL(TE1)={il_te1:.4f} dB"
    )

    if (i + 1) % 15 == 0 or i == 0:
        sim_viz = get_sim(params, beta=beta, input_mode=0)
        _, ax = plt.subplots(figsize=(5, 3))
        sim_viz.plot_eps(z=0.01, ax=ax, monitor_alpha=0.0, source_alpha=0.0)
        ax.set_title(f"Step {i+1} [{stage}]")
        plt.show()

_fab_weight[0] = 1.0
params_standard_to = param_history_p1[-1]
print("\nPhase 1 complete.")
======================================================================
PHASE 1: Staged TO (unconstrained → ramp → full fab)
======================================================================












  Step   1/75, β=8.0, fab_w=0.00 [unconstrained]  obj=0.095115  IL(TE0)=8.6810 dB  IL(TE1)=12.6168 dB
































  Step   2/75, β=8.3, fab_w=0.00 [unconstrained]  obj=0.425841  IL(TE0)=1.8478 dB  IL(TE1)=7.0287 dB


















  Step   3/75, β=8.6, fab_w=0.00 [unconstrained]  obj=0.471936  IL(TE0)=1.3040 dB  IL(TE1)=6.9200 dB


















  Step   4/75, β=9.0, fab_w=0.00 [unconstrained]  obj=0.511760  IL(TE0)=1.3737 dB  IL(TE1)=5.3064 dB


















  Step   5/75, β=9.3, fab_w=0.00 [unconstrained]  obj=0.676725  IL(TE0)=0.5659 dB  IL(TE1)=3.2274 dB


















  Step   6/75, β=9.6, fab_w=0.00 [unconstrained]  obj=0.741746  IL(TE0)=0.8295 dB  IL(TE1)=1.8220 dB


















  Step   7/75, β=9.9, fab_w=0.00 [unconstrained]  obj=0.780133  IL(TE0)=0.6412 dB  IL(TE1)=1.5644 dB


















  Step   8/75, β=10.3, fab_w=0.00 [unconstrained]  obj=0.817891  IL(TE0)=0.6547 dB  IL(TE1)=1.1029 dB


















  Step   9/75, β=10.6, fab_w=0.00 [unconstrained]  obj=0.875817  IL(TE0)=0.4423 dB  IL(TE1)=0.7137 dB


















  Step  10/75, β=10.9, fab_w=0.00 [unconstrained]  obj=0.897636  IL(TE0)=0.3046 dB  IL(TE1)=0.6399 dB


















  Step  11/75, β=11.2, fab_w=0.00 [unconstrained]  obj=0.883930  IL(TE0)=0.4420 dB  IL(TE1)=0.6317 dB


















  Step  12/75, β=11.6, fab_w=0.00 [unconstrained]  obj=0.921600  IL(TE0)=0.2621 dB  IL(TE1)=0.4490 dB


















  Step  13/75, β=11.9, fab_w=0.00 [unconstrained]  obj=0.936383  IL(TE0)=0.1622 dB  IL(TE1)=0.4123 dB


















  Step  14/75, β=12.2, fab_w=0.00 [unconstrained]  obj=0.941380  IL(TE0)=0.2236 dB  IL(TE1)=0.3015 dB


















  Step  15/75, β=12.5, fab_w=0.00 [unconstrained]  obj=0.957070  IL(TE0)=0.1462 dB  IL(TE1)=0.2354 dB
































  Step  16/75, β=12.9, fab_w=0.00 [unconstrained]  obj=0.942088  IL(TE0)=0.2112 dB  IL(TE1)=0.3075 dB


















  Step  17/75, β=13.2, fab_w=0.00 [unconstrained]  obj=0.957995  IL(TE0)=0.1270 dB  IL(TE1)=0.2465 dB


















  Step  18/75, β=13.5, fab_w=0.00 [unconstrained]  obj=0.960078  IL(TE0)=0.1152 dB  IL(TE1)=0.2395 dB


















  Step  19/75, β=13.8, fab_w=0.00 [unconstrained]  obj=0.963912  IL(TE0)=0.1296 dB  IL(TE1)=0.1899 dB


















  Step  20/75, β=14.2, fab_w=0.00 [unconstrained]  obj=0.971937  IL(TE0)=0.0840 dB  IL(TE1)=0.1636 dB


















  Step  21/75, β=14.5, fab_w=0.00 [unconstrained]  obj=0.970265  IL(TE0)=0.1265 dB  IL(TE1)=0.1357 dB


















  Step  22/75, β=14.8, fab_w=0.00 [unconstrained]  obj=0.978262  IL(TE0)=0.0725 dB  IL(TE1)=0.1185 dB


















  Step  23/75, β=15.1, fab_w=0.00 [unconstrained]  obj=0.975426  IL(TE0)=0.0913 dB  IL(TE1)=0.1248 dB


















  Step  24/75, β=15.5, fab_w=0.00 [unconstrained]  obj=0.981541  IL(TE0)=0.0661 dB  IL(TE1)=0.0958 dB


















  Step  25/75, β=15.8, fab_w=0.00 [unconstrained]  obj=0.983202  IL(TE0)=0.0680 dB  IL(TE1)=0.0791 dB


















  Step  26/75, β=16.1, fab_w=0.00 [ramping]  obj=0.981251  IL(TE0)=0.0728 dB  IL(TE1)=0.0916 dB


















  Step  27/75, β=16.4, fab_w=0.03 [ramping]  obj=0.968142  IL(TE0)=0.0596 dB  IL(TE1)=0.0894 dB


















  Step  28/75, β=16.8, fab_w=0.05 [ramping]  obj=0.954094  IL(TE0)=0.0714 dB  IL(TE1)=0.0761 dB


















  Step  29/75, β=17.1, fab_w=0.07 [ramping]  obj=0.946372  IL(TE0)=0.0472 dB  IL(TE1)=0.0486 dB


















  Step  30/75, β=17.4, fab_w=0.10 [ramping]  obj=0.931373  IL(TE0)=0.0628 dB  IL(TE1)=0.0513 dB
































  Step  31/75, β=17.7, fab_w=0.12 [ramping]  obj=0.922108  IL(TE0)=0.0399 dB  IL(TE1)=0.0482 dB


















  Step  32/75, β=18.1, fab_w=0.15 [ramping]  obj=0.910643  IL(TE0)=0.0531 dB  IL(TE1)=0.0357 dB


















  Step  33/75, β=18.4, fab_w=0.17 [ramping]  obj=0.901825  IL(TE0)=0.0367 dB  IL(TE1)=0.0379 dB


















  Step  34/75, β=18.7, fab_w=0.20 [ramping]  obj=0.891365  IL(TE0)=0.0437 dB  IL(TE1)=0.0394 dB


















  Step  35/75, β=19.0, fab_w=0.23 [ramping]  obj=0.884377  IL(TE0)=0.0345 dB  IL(TE1)=0.0360 dB


















  Step  36/75, β=19.4, fab_w=0.25 [ramping]  obj=0.877490  IL(TE0)=0.0334 dB  IL(TE1)=0.0329 dB


















  Step  37/75, β=19.7, fab_w=0.28 [ramping]  obj=0.872901  IL(TE0)=0.0298 dB  IL(TE1)=0.0207 dB


















  Step  38/75, β=20.0, fab_w=0.30 [ramping]  obj=0.867207  IL(TE0)=0.0284 dB  IL(TE1)=0.0255 dB


















  Step  39/75, β=20.3, fab_w=0.33 [ramping]  obj=0.863982  IL(TE0)=0.0260 dB  IL(TE1)=0.0193 dB


















  Step  40/75, β=20.6, fab_w=0.35 [ramping]  obj=0.860994  IL(TE0)=0.0255 dB  IL(TE1)=0.0168 dB


















  Step  41/75, β=21.0, fab_w=0.38 [ramping]  obj=0.858790  IL(TE0)=0.0238 dB  IL(TE1)=0.0174 dB


















  Step  42/75, β=21.3, fab_w=0.40 [ramping]  obj=0.857432  IL(TE0)=0.0296 dB  IL(TE1)=0.0125 dB


















  Step  43/75, β=21.6, fab_w=0.42 [ramping]  obj=0.857701  IL(TE0)=0.0250 dB  IL(TE1)=0.0100 dB


















  Step  44/75, β=21.9, fab_w=0.45 [ramping]  obj=0.857920  IL(TE0)=0.0222 dB  IL(TE1)=0.0120 dB


















  Step  45/75, β=22.3, fab_w=0.47 [ramping]  obj=0.858758  IL(TE0)=0.0218 dB  IL(TE1)=0.0148 dB
































  Step  46/75, β=22.6, fab_w=0.50 [ramping]  obj=0.861020  IL(TE0)=0.0201 dB  IL(TE1)=0.0131 dB


















  Step  47/75, β=22.9, fab_w=0.53 [ramping]  obj=0.863491  IL(TE0)=0.0212 dB  IL(TE1)=0.0121 dB


















  Step  48/75, β=23.2, fab_w=0.55 [ramping]  obj=0.867596  IL(TE0)=0.0178 dB  IL(TE1)=0.0094 dB


















  Step  49/75, β=23.6, fab_w=0.57 [ramping]  obj=0.872126  IL(TE0)=0.0196 dB  IL(TE1)=0.0087 dB


















  Step  50/75, β=23.9, fab_w=0.60 [ramping]  obj=0.877990  IL(TE0)=0.0196 dB  IL(TE1)=0.0078 dB


















  Step  51/75, β=24.2, fab_w=0.62 [ramping]  obj=0.883549  IL(TE0)=0.0231 dB  IL(TE1)=0.0071 dB


















  Step  52/75, β=24.5, fab_w=0.65 [ramping]  obj=0.889408  IL(TE0)=0.0220 dB  IL(TE1)=0.0055 dB


















  Step  53/75, β=24.9, fab_w=0.68 [ramping]  obj=0.894928  IL(TE0)=0.0231 dB  IL(TE1)=0.0052 dB


















  Step  54/75, β=25.2, fab_w=0.70 [ramping]  obj=0.900780  IL(TE0)=0.0243 dB  IL(TE1)=0.0047 dB


















  Step  55/75, β=25.5, fab_w=0.72 [ramping]  obj=0.905956  IL(TE0)=0.0264 dB  IL(TE1)=0.0059 dB


















  Step  56/75, β=25.8, fab_w=0.75 [ramping]  obj=0.910700  IL(TE0)=0.0262 dB  IL(TE1)=0.0064 dB


















  Step  57/75, β=26.2, fab_w=0.78 [ramping]  obj=0.914859  IL(TE0)=0.0263 dB  IL(TE1)=0.0050 dB


















  Step  58/75, β=26.5, fab_w=0.80 [ramping]  obj=0.917980  IL(TE0)=0.0270 dB  IL(TE1)=0.0054 dB


















  Step  59/75, β=26.8, fab_w=0.82 [ramping]  obj=0.920749  IL(TE0)=0.0272 dB  IL(TE1)=0.0065 dB


















  Step  60/75, β=27.1, fab_w=0.85 [ramping]  obj=0.922851  IL(TE0)=0.0291 dB  IL(TE1)=0.0081 dB
































  Step  61/75, β=27.5, fab_w=0.88 [ramping]  obj=0.924669  IL(TE0)=0.0283 dB  IL(TE1)=0.0075 dB


















  Step  62/75, β=27.8, fab_w=0.90 [ramping]  obj=0.925981  IL(TE0)=0.0274 dB  IL(TE1)=0.0069 dB


















  Step  63/75, β=28.1, fab_w=0.93 [ramping]  obj=0.926952  IL(TE0)=0.0283 dB  IL(TE1)=0.0078 dB


















  Step  64/75, β=28.4, fab_w=0.95 [ramping]  obj=0.927686  IL(TE0)=0.0291 dB  IL(TE1)=0.0094 dB


















  Step  65/75, β=28.8, fab_w=0.97 [ramping]  obj=0.928128  IL(TE0)=0.0298 dB  IL(TE1)=0.0113 dB


















  Step  66/75, β=29.1, fab_w=1.00 [full-fab]  obj=0.928452  IL(TE0)=0.0293 dB  IL(TE1)=0.0121 dB


















  Step  67/75, β=29.4, fab_w=1.00 [full-fab]  obj=0.930192  IL(TE0)=0.0282 dB  IL(TE1)=0.0120 dB


















  Step  68/75, β=29.7, fab_w=1.00 [full-fab]  obj=0.931774  IL(TE0)=0.0287 dB  IL(TE1)=0.0124 dB


















  Step  69/75, β=30.1, fab_w=1.00 [full-fab]  obj=0.933152  IL(TE0)=0.0298 dB  IL(TE1)=0.0135 dB


















  Step  70/75, β=30.4, fab_w=1.00 [full-fab]  obj=0.934433  IL(TE0)=0.0294 dB  IL(TE1)=0.0140 dB


















  Step  71/75, β=30.7, fab_w=1.00 [full-fab]  obj=0.935582  IL(TE0)=0.0289 dB  IL(TE1)=0.0148 dB


















  Step  72/75, β=31.0, fab_w=1.00 [full-fab]  obj=0.936588  IL(TE0)=0.0286 dB  IL(TE1)=0.0157 dB


















  Step  73/75, β=31.4, fab_w=1.00 [full-fab]  obj=0.937569  IL(TE0)=0.0293 dB  IL(TE1)=0.0164 dB


















  Step  74/75, β=31.7, fab_w=1.00 [full-fab]  obj=0.938437  IL(TE0)=0.0306 dB  IL(TE1)=0.0169 dB


















  Step  75/75, β=32.0, fab_w=1.00 [full-fab]  obj=0.939263  IL(TE0)=0.0299 dB  IL(TE1)=0.0165 dB





















Phase 1 complete.










# Phase 1 IL convergence (per-mode)
fig, ax = plt.subplots(figsize=(9, 4))
ax.plot(
    range(1, len(il_te0_history_p1) + 1), il_te0_history_p1, "b-", label="IL(TE0 → bot)"
)
ax.plot(
    range(1, len(il_te1_history_p1) + 1), il_te1_history_p1, "r-", label="IL(TE1 → top)"
)
ax.plot(
    range(1, len(il_history_p1) + 1),
    il_history_p1,
    "k--",
    linewidth=1.5,
    alpha=0.7,
    label="Average IL",
)
ax.axvline(x=unconstrained_end, color="gray", linestyle=":", alpha=0.6)
ax.axvline(x=ramp_end, color="gray", linestyle=":", alpha=0.6)
ax.text(
    unconstrained_end / 2,
    ax.get_ylim()[1] * 0.95,
    "unconstrained",
    ha="center",
    fontsize=9,
)
ax.text(
    (unconstrained_end + ramp_end) / 2,
    ax.get_ylim()[1] * 0.95,
    "ramping",
    ha="center",
    fontsize=9,
)
ax.text(
    (ramp_end + num_steps_phase1) / 2,
    ax.get_ylim()[1] * 0.95,
    "full-fab",
    ha="center",
    fontsize=9,
)
ax.set_xlabel("Iteration")
ax.set_ylabel("Insertion Loss (dB)")
ax.set_title("Phase 1 Convergence — Per-Mode and Average IL")
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

























7. Phase 2: Two Parallel Branches From the Iter-135 Design


Both branches start from the same Phase-1 final params and run the same number of additional Adam steps.


    

  • 
Branch A (standard TO) simply continues at fixed β = 32.
    • 

  • 
Branch B (contour-constrained) re-softens the projection and ramps β from 12 → 32 while the contour penalty is active. This is the key fix: at the terminal β = 32 the tanh projection is saturated (∂density/∂params ≈ 0), so the contour penalty has no gradient leverage no matter how large its weight. Applying it while β is still ramping keeps parasitic features erodible, so the penalty can actually drive them to zero.
    • 



The contour mask is recomputed every mask_update_freq steps from a forward simulation of the TE0-in sub-problem (the dominant routing pathway); any connected feature carrying negligible field there is marked for erosion.






T_P = 0.05  # mark contours holding less than 5% of total field power
contour_weight = 0.75
mask_update_freq = 10
pixel_area = pixel_size**2

# Field monitor over optimization frequencies, only used for contour analysis
field_monitor_contour = td.FieldMonitor(
    center=(0, 0, 0),
    size=(td.inf, td.inf, 0),
    freqs=freqs_opt,
    name="field",
)


def compute_contour_penalty_mask(
    params: np.ndarray, beta: float, threshold: float
) -> np.ndarray:
    """Compute the indicator mask of voxels in low-field contours.

    Uses the TE0-in forward simulation: this is the dominant routing pathway
    and any feature with negligible field there can be safely eroded."""
    density = get_density(params, beta=beta).reshape(num_pixels_x, num_pixels_y)
    binary_design = (density > 0.5).astype(int)
    labeled_contours, num_contours_local = scipy.ndimage.label(binary_design)
    n_map = np.where(binary_design, n_si, n_sio2)

    sim_field = get_sim(params, beta=beta, input_mode=0)
    sim_field = sim_field.updated_copy(
        monitors=list(sim_field.monitors) + [field_monitor_contour]
    )
    sim_data_field = web.run(
        sim_field,
        folder_name="contour constraint demux",
        task_name="contour_analysis",
        verbose=False,
    )
    field_data = sim_data_field["field"]
    E_sq_per_freq = (
        np.abs(field_data.Ex.values) ** 2
        + np.abs(field_data.Ey.values) ** 2
        + np.abs(field_data.Ez.values) ** 2
    )
    E_sq = np.mean(E_sq_per_freq, axis=-1)

    x_field = field_data.Ex.coords["x"].values
    y_field = field_data.Ex.coords["y"].values
    x_design_grid = np.linspace(
        -design_size_x / 2 + pixel_size / 2,
        design_size_x / 2 - pixel_size / 2,
        num_pixels_x,
    )
    y_design_grid = np.linspace(
        -design_size_y / 2 + pixel_size / 2,
        design_size_y / 2 - pixel_size / 2,
        num_pixels_y,
    )

    interp = RegularGridInterpolator(
        (x_field, y_field),
        E_sq[:, :, 0],
        method="linear",
        bounds_error=False,
        fill_value=0.0,
    )
    xx_des, yy_des = np.meshgrid(x_design_grid, y_design_grid, indexing="ij")
    E_sq_design_local = interp(
        np.stack([xx_des.ravel(), yy_des.ravel()], axis=-1)
    ).reshape(num_pixels_x, num_pixels_y)

    P_tot_local = np.sum(n_map * E_sq_design_local)
    mask_out = np.zeros((num_pixels_x, num_pixels_y, 1))
    n_marked = 0
    for m in range(1, num_contours_local + 1):
        mask_m = labeled_contours == m
        P_m = np.sum(np.where(mask_m, n_map, 0) * E_sq_design_local)
        if P_tot_local > 0 and P_m / P_tot_local < threshold:
            mask_out += mask_m.reshape(num_pixels_x, num_pixels_y, 1).astype(float)
            n_marked += 1
    print(
        f"    [contour mask] {num_contours_local} contours, {n_marked} marked (P/P_tot < {threshold:.2f})"
    )
    return mask_out


_contour_mask_holder = [np.zeros((num_pixels_x, num_pixels_y, 1))]


def objective_with_contour(params: np.ndarray, beta: float) -> float:
    transmission = get_transmission(params, beta=beta)
    fab_penalty = get_penalty(params, beta=beta)
    density = get_density(params, beta=beta)
    contour_penalty = np.sum(_contour_mask_holder[0] * density) * pixel_area
    return (
        transmission - _fab_weight[0] * fab_penalty - contour_weight * contour_penalty
    )











num_steps_phase2 = 25
learning_rate_p2 = 0.02
beta_p2_start = 12  # re-soften so the contour penalty can actually erode features
_fab_weight[0] = 1.0  # full fab strength for all later phases


def get_beta_p2(step_num: int, num_steps: int, use_contour: bool) -> float:
    """Contour branch ramps β so features stay erodible while the penalty acts;
    the standard baseline just continues at full β."""
    if not use_contour:
        return beta_max
    return beta_p2_start + (beta_max - beta_p2_start) * step_num / max(num_steps - 1, 1)


def run_phase2(params_init: np.ndarray, use_contour: bool, label: str):
    params = params_init.copy()
    optimizer = AdamOptimizer.construct(learning_rate=learning_rate_p2)
    opt_state = optimizer.initial_state(params)

    beta0 = get_beta_p2(0, num_steps_phase2, use_contour)
    if use_contour:
        _contour_mask_holder[0] = compute_contour_penalty_mask(
            params, beta0, threshold=T_P
        )
        val_grad_fn_local = autograd.value_and_grad(objective_with_contour)
    else:
        val_grad_fn_local = autograd.value_and_grad(objective)

    history_params = [params.copy()]
    history_il_te0 = []
    history_il_te1 = []

    print("=" * 70)
    print(f"PHASE 2 — {label}")
    print("=" * 70)

    for i in range(num_steps_phase2):
        beta = get_beta_p2(i, num_steps_phase2, use_contour)
        if use_contour and i > 0 and i % mask_update_freq == 0:
            _contour_mask_holder[0] = compute_contour_penalty_mask(
                params, beta, threshold=T_P
            )

        value, gradient = val_grad_fn_local(params, beta=beta)
        params, opt_state = optimizer.update(
            parameters=params, gradient=-gradient, state=opt_state
        )
        params = np.clip(params, 0, 1)

        il_te0 = -10 * np.log10(max(_last_T_te0[0], 1e-12))
        il_te1 = -10 * np.log10(max(_last_T_te1[0], 1e-12))
        history_il_te0.append(il_te0)
        history_il_te1.append(il_te1)
        history_params.append(params.copy())
        print(
            f"  Step {i+1:3d}/{num_steps_phase2}  β={beta:.1f}  obj={value:.6f}  IL(TE0)={il_te0:.4f} dB  IL(TE1)={il_te1:.4f} dB"
        )

        if (i + 1) % 20 == 0 or i == 0:
            sim_viz = get_sim(params, beta=beta, input_mode=0)
            _, ax = plt.subplots(figsize=(5, 3))
            sim_viz.plot_eps(z=0.01, ax=ax, monitor_alpha=0.0, source_alpha=0.0)
            ax.set_title(f"{label}, Step {i+1}")
            plt.show()

    print(f"\n{label} complete.")
    return history_params, history_il_te0, history_il_te1


# Branch A: continue without contour constraint
param_history_p2_std, il_te0_p2_std, il_te1_p2_std = run_phase2(
    params_standard_to, use_contour=False, label="Standard-TO Continued"
)

# Branch B: with the contour constraint, applied while β ramps 12 -> 32
param_history_p2_cc, il_te0_p2_cc, il_te1_p2_cc = run_phase2(
    params_standard_to, use_contour=True, label="Contour-Constrained"
)







======================================================================
PHASE 2 — Standard-TO Continued
======================================================================
/var/folders/qn/syhrzy8n7930sgqvxv2x65s40000gn/T/ipykernel_62177/3064855979.py:17: PydanticDeprecatedSince20: The `construct` method is deprecated; use `model_construct` instead. Deprecated in Pydantic V2.0 to be removed in V3.0. See Pydantic V2 Migration Guide at https://errors.pydantic.dev/2.10/migration/
  optimizer = AdamOptimizer.construct(learning_rate=learning_rate_p2)


















  Step   1/25  β=32.0  obj=0.939977  IL(TE0)=0.0301 dB  IL(TE1)=0.0168 dB
































  Step   2/25  β=32.0  obj=0.770512  IL(TE0)=0.6002 dB  IL(TE1)=1.0682 dB


















  Step   3/25  β=32.0  obj=0.909626  IL(TE0)=0.1182 dB  IL(TE1)=0.2605 dB


















  Step   4/25  β=32.0  obj=0.841428  IL(TE0)=0.3869 dB  IL(TE1)=0.6714 dB


















  Step   5/25  β=32.0  obj=0.889329  IL(TE0)=0.2401 dB  IL(TE1)=0.3681 dB


















  Step   6/25  β=32.0  obj=0.923377  IL(TE0)=0.0916 dB  IL(TE1)=0.1785 dB


















  Step   7/25  β=32.0  obj=0.907351  IL(TE0)=0.0928 dB  IL(TE1)=0.2972 dB


















  Step   8/25  β=32.0  obj=0.898345  IL(TE0)=0.1186 dB  IL(TE1)=0.3481 dB


















  Step   9/25  β=32.0  obj=0.915820  IL(TE0)=0.0944 dB  IL(TE1)=0.2226 dB


















  Step  10/25  β=32.0  obj=0.934899  IL(TE0)=0.0606 dB  IL(TE1)=0.1065 dB


















  Step  11/25  β=32.0  obj=0.933037  IL(TE0)=0.0703 dB  IL(TE1)=0.1329 dB


















  Step  12/25  β=32.0  obj=0.922296  IL(TE0)=0.1026 dB  IL(TE1)=0.2061 dB


















  Step  13/25  β=32.0  obj=0.922429  IL(TE0)=0.1079 dB  IL(TE1)=0.1965 dB


















  Step  14/25  β=32.0  obj=0.932617  IL(TE0)=0.0797 dB  IL(TE1)=0.1231 dB


















  Step  15/25  β=32.0  obj=0.941455  IL(TE0)=0.0475 dB  IL(TE1)=0.0673 dB


















  Step  16/25  β=32.0  obj=0.941677  IL(TE0)=0.0407 dB  IL(TE1)=0.0657 dB


















  Step  17/25  β=32.0  obj=0.936945  IL(TE0)=0.0543 dB  IL(TE1)=0.0920 dB


















  Step  18/25  β=32.0  obj=0.935165  IL(TE0)=0.0617 dB  IL(TE1)=0.1036 dB


















  Step  19/25  β=32.0  obj=0.938844  IL(TE0)=0.0528 dB  IL(TE1)=0.0878 dB


















  Step  20/25  β=32.0  obj=0.943887  IL(TE0)=0.0416 dB  IL(TE1)=0.0647 dB
































  Step  21/25  β=32.0  obj=0.945723  IL(TE0)=0.0435 dB  IL(TE1)=0.0555 dB


















  Step  22/25  β=32.0  obj=0.944118  IL(TE0)=0.0565 dB  IL(TE1)=0.0596 dB


















  Step  23/25  β=32.0  obj=0.943374  IL(TE0)=0.0645 dB  IL(TE1)=0.0619 dB


















  Step  24/25  β=32.0  obj=0.944666  IL(TE0)=0.0582 dB  IL(TE1)=0.0551 dB


















  Step  25/25  β=32.0  obj=0.946774  IL(TE0)=0.0438 dB  IL(TE1)=0.0460 dB

Standard-TO Continued complete.
/var/folders/qn/syhrzy8n7930sgqvxv2x65s40000gn/T/ipykernel_62177/3064855979.py:17: PydanticDeprecatedSince20: The `construct` method is deprecated; use `model_construct` instead. Deprecated in Pydantic V2.0 to be removed in V3.0. See Pydantic V2 Migration Guide at https://errors.pydantic.dev/2.10/migration/
  optimizer = AdamOptimizer.construct(learning_rate=learning_rate_p2)
    [contour mask] 4 contours, 3 marked (P/P_tot < 0.05)
======================================================================
PHASE 2 — Contour-Constrained
======================================================================


















  Step   1/25  β=12.0  obj=-4.306218  IL(TE0)=0.3976 dB  IL(TE1)=0.5276 dB
































  Step   2/25  β=12.8  obj=-4.231834  IL(TE0)=0.6491 dB  IL(TE1)=1.0511 dB


















  Step   3/25  β=13.7  obj=-3.971378  IL(TE0)=0.2631 dB  IL(TE1)=0.5384 dB


















  Step   4/25  β=14.5  obj=-3.738763  IL(TE0)=0.1895 dB  IL(TE1)=0.2860 dB


















  Step   5/25  β=15.3  obj=-3.470009  IL(TE0)=0.2446 dB  IL(TE1)=0.3031 dB


















  Step   6/25  β=16.2  obj=-3.069188  IL(TE0)=0.1997 dB  IL(TE1)=0.2714 dB


















  Step   7/25  β=17.0  obj=-2.542798  IL(TE0)=0.1284 dB  IL(TE1)=0.1565 dB


















  Step   8/25  β=17.8  obj=-1.936563  IL(TE0)=0.1359 dB  IL(TE1)=0.1052 dB


















  Step   9/25  β=18.7  obj=-1.289786  IL(TE0)=0.1733 dB  IL(TE1)=0.1390 dB


















  Step  10/25  β=19.5  obj=-0.688261  IL(TE0)=0.1438 dB  IL(TE1)=0.1449 dB
    [contour mask] 29 contours, 28 marked (P/P_tot < 0.05)


















  Step  11/25  β=20.3  obj=0.245071  IL(TE0)=0.0862 dB  IL(TE1)=0.1142 dB


















  Step  12/25  β=21.2  obj=0.393050  IL(TE0)=0.0655 dB  IL(TE1)=0.0935 dB


















  Step  13/25  β=22.0  obj=0.547862  IL(TE0)=0.0804 dB  IL(TE1)=0.0889 dB


















  Step  14/25  β=22.8  obj=0.690882  IL(TE0)=0.0854 dB  IL(TE1)=0.0737 dB


















  Step  15/25  β=23.7  obj=0.803066  IL(TE0)=0.0753 dB  IL(TE1)=0.0596 dB


















  Step  16/25  β=24.5  obj=0.870381  IL(TE0)=0.0679 dB  IL(TE1)=0.0617 dB


















  Step  17/25  β=25.3  obj=0.903139  IL(TE0)=0.0640 dB  IL(TE1)=0.0614 dB


















  Step  18/25  β=26.2  obj=0.921612  IL(TE0)=0.0560 dB  IL(TE1)=0.0484 dB


















  Step  19/25  β=27.0  obj=0.930805  IL(TE0)=0.0496 dB  IL(TE1)=0.0468 dB


















  Step  20/25  β=27.8  obj=0.933617  IL(TE0)=0.0468 dB  IL(TE1)=0.0571 dB




















    [contour mask] 1 contours, 0 marked (P/P_tot < 0.05)


















  Step  21/25  β=28.7  obj=0.935255  IL(TE0)=0.0453 dB  IL(TE1)=0.0556 dB


















  Step  22/25  β=29.5  obj=0.934245  IL(TE0)=0.0529 dB  IL(TE1)=0.0493 dB


















  Step  23/25  β=30.3  obj=0.933565  IL(TE0)=0.0597 dB  IL(TE1)=0.0449 dB


















  Step  24/25  β=31.2  obj=0.936553  IL(TE0)=0.0503 dB  IL(TE1)=0.0340 dB


















  Step  25/25  β=32.0  obj=0.940257  IL(TE0)=0.0391 dB  IL(TE1)=0.0273 dB

Contour-Constrained complete.










# Per-mode Phase-1 + Phase-2 convergence comparison
fig, axes = plt.subplots(2, 1, figsize=(10, 7), sharex=True)

n_p1 = len(il_te0_history_p1)
n_p2 = len(il_te0_p2_std)
x_p1 = np.arange(1, n_p1 + 1)
x_p2 = np.arange(n_p1 + 1, n_p1 + n_p2 + 1)

axes[0].plot(x_p1, il_te0_history_p1, "k-", linewidth=1.5, label="Phase 1 (shared)")
axes[0].plot(
    x_p2, il_te0_p2_std, "b-", linewidth=1.5, label="Phase 2: Std-TO continued"
)
axes[0].plot(
    x_p2, il_te0_p2_cc, "r-", linewidth=1.5, label="Phase 2: Contour-constrained"
)
axes[0].axvline(x=n_p1 + 0.5, color="gray", linestyle="--", alpha=0.6)
axes[0].set_ylabel("IL (TE0 → bot, dB)")
axes[0].set_title("TE0 routing — Phase 1 + Phase 2")
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(x_p1, il_te1_history_p1, "k-", linewidth=1.5, label="Phase 1 (shared)")
axes[1].plot(
    x_p2, il_te1_p2_std, "b-", linewidth=1.5, label="Phase 2: Std-TO continued"
)
axes[1].plot(
    x_p2, il_te1_p2_cc, "r-", linewidth=1.5, label="Phase 2: Contour-constrained"
)
axes[1].axvline(x=n_p1 + 0.5, color="gray", linestyle="--", alpha=0.6)
axes[1].set_xlabel("Iteration")
axes[1].set_ylabel("IL (TE1 → top, dB)")
axes[1].set_title("TE1 mode-conversion + routing — Phase 1 + Phase 2")
axes[1].legend()
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

























8. Phase 3: Seeded-TO Finishing


Both branches receive a fab-only seeded-TO polishing pass. Critically, the contour penalty is off during this stage — it has done its job and should not keep eroding material, which would damage the main contour. This matches the paper’s description that seeded TO smooths rough edges and enforces foundry DRC.






num_steps_seeded = 10
learning_rate_seeded = 0.005


def run_seeded_to(params_init: np.ndarray, label: str):
    """Polishing pass: small LR, fab-only objective, soft re-seed."""
    density_seed = get_density(params_init, beta=beta_max).reshape(
        num_pixels_x, num_pixels_y, 1
    )
    binary_seed = (density_seed > 0.5).astype(float)
    params = 0.05 + 0.9 * binary_seed

    optimizer = AdamOptimizer.construct(learning_rate=learning_rate_seeded)
    opt_state = optimizer.initial_state(params)

    val_grad_fn_local = autograd.value_and_grad(objective)

    history_il_te0 = []
    history_il_te1 = []
    print("=" * 70)
    print(f"PHASE 3: Seeded TO — {label} (fab-only)")
    print("=" * 70)

    for i in range(num_steps_seeded):
        value, gradient = val_grad_fn_local(params, beta=beta_max)
        params, opt_state = optimizer.update(
            parameters=params, gradient=-gradient, state=opt_state
        )
        params = np.clip(params, 0, 1)

        il_te0 = -10 * np.log10(max(_last_T_te0[0], 1e-12))
        il_te1 = -10 * np.log10(max(_last_T_te1[0], 1e-12))
        history_il_te0.append(il_te0)
        history_il_te1.append(il_te1)
        print(
            f"  Step {i+1:3d}/{num_steps_seeded}  obj={value:.6f}  IL(TE0)={il_te0:.4f} dB  IL(TE1)={il_te1:.4f} dB"
        )

    return params, history_il_te0, history_il_te1


_fab_weight[0] = 1.0
params_standard_seeded, il_te0_seed_std, il_te1_seed_std = run_seeded_to(
    param_history_p2_std[-1], label="Standard-TO Continued"
)
params_contour_seeded, il_te0_seed_cc, il_te1_seed_cc = run_seeded_to(
    param_history_p2_cc[-1], label="Contour-Constrained"
)







======================================================================
PHASE 3: Seeded TO — Standard-TO Continued (fab-only)
======================================================================
/var/folders/qn/syhrzy8n7930sgqvxv2x65s40000gn/T/ipykernel_62177/3601026129.py:13: PydanticDeprecatedSince20: The `construct` method is deprecated; use `model_construct` instead. Deprecated in Pydantic V2.0 to be removed in V3.0. See Pydantic V2 Migration Guide at https://errors.pydantic.dev/2.10/migration/
  optimizer = AdamOptimizer.construct(learning_rate=learning_rate_seeded)


















  Step   1/10  obj=0.614772  IL(TE0)=1.7591 dB  IL(TE1)=1.7964 dB


















  Step   2/10  obj=0.634819  IL(TE0)=1.6368 dB  IL(TE1)=1.6594 dB


















  Step   3/10  obj=0.654505  IL(TE0)=1.5192 dB  IL(TE1)=1.5294 dB


















  Step   4/10  obj=0.673859  IL(TE0)=1.4056 dB  IL(TE1)=1.4060 dB


















  Step   5/10  obj=0.692915  IL(TE0)=1.2952 dB  IL(TE1)=1.2885 dB


















  Step   6/10  obj=0.711697  IL(TE0)=1.1877 dB  IL(TE1)=1.1768 dB


















  Step   7/10  obj=0.730217  IL(TE0)=1.0824 dB  IL(TE1)=1.0704 dB


















  Step   8/10  obj=0.748477  IL(TE0)=0.9790 dB  IL(TE1)=0.9694 dB


















  Step   9/10  obj=0.766456  IL(TE0)=0.8770 dB  IL(TE1)=0.8740 dB


















  Step  10/10  obj=0.784059  IL(TE0)=0.7763 dB  IL(TE1)=0.7850 dB
======================================================================
PHASE 3: Seeded TO — Contour-Constrained (fab-only)
======================================================================
/var/folders/qn/syhrzy8n7930sgqvxv2x65s40000gn/T/ipykernel_62177/3601026129.py:13: PydanticDeprecatedSince20: The `construct` method is deprecated; use `model_construct` instead. Deprecated in Pydantic V2.0 to be removed in V3.0. See Pydantic V2 Migration Guide at https://errors.pydantic.dev/2.10/migration/
  optimizer = AdamOptimizer.construct(learning_rate=learning_rate_seeded)


















  Step   1/10  obj=0.661275  IL(TE0)=1.6087 dB  IL(TE1)=1.4612 dB


















  Step   2/10  obj=0.680328  IL(TE0)=1.4694 dB  IL(TE1)=1.3633 dB


















  Step   3/10  obj=0.699242  IL(TE0)=1.3345 dB  IL(TE1)=1.2693 dB


















  Step   4/10  obj=0.718016  IL(TE0)=1.2037 dB  IL(TE1)=1.1789 dB


















  Step   5/10  obj=0.736598  IL(TE0)=1.0772 dB  IL(TE1)=1.0920 dB


















  Step   6/10  obj=0.754893  IL(TE0)=0.9558 dB  IL(TE1)=1.0083 dB


















  Step   7/10  obj=0.772767  IL(TE0)=0.8402 dB  IL(TE1)=0.9279 dB


















  Step   8/10  obj=0.790041  IL(TE0)=0.7314 dB  IL(TE1)=0.8508 dB


















  Step   9/10  obj=0.806475  IL(TE0)=0.6306 dB  IL(TE1)=0.7774 dB


















  Step  10/10  obj=0.821767  IL(TE0)=0.5390 dB  IL(TE1)=0.7084 dB










9. Broadband Comparison (1260–1360 nm)


Replicates Fig. 4c of the paper: per-mode insertion loss spectra for both the standard-TO and contour-constrained final devices.






def run_broadband(params_final: np.ndarray, label: str):
    """Run TE0-in and TE1-in broadband sims; return (T_te0, T_te1) over freqs_analysis."""
    sim_te0 = get_sim(
        params_final, beta=beta_max, input_mode=0, broadband=True, with_field=True
    )
    sim_te1 = get_sim(
        params_final, beta=beta_max, input_mode=1, broadband=True, with_field=True
    )

    batch = web.run_async(
        {f"{label}_te0": sim_te0, f"{label}_te1": sim_te1},
        folder_name="contour constraint demux",
        verbose=True
    )

    amp_te0_to_bot = (
        batch[f"{label}_te0"]["mode_bot"]
        .amps.sel(mode_index=0, direction="+")
        .values.squeeze()
    )
    amp_te1_to_top = (
        batch[f"{label}_te1"]["mode_top"]
        .amps.sel(mode_index=0, direction="+")
        .values.squeeze()
    )

    return np.abs(amp_te0_to_bot) ** 2, np.abs(amp_te1_to_top) ** 2, batch


T_std_te0, T_std_te1, batch_std = run_broadband(params_standard_seeded, "std")
T_cc_te0, T_cc_te1, batch_cc = run_broadband(params_contour_seeded, "cc")

IL_std_te0 = -10 * np.log10(np.maximum(T_std_te0, 1e-12))
IL_std_te1 = -10 * np.log10(np.maximum(T_std_te1, 1e-12))
IL_cc_te0 = -10 * np.log10(np.maximum(T_cc_te0, 1e-12))
IL_cc_te1 = -10 * np.log10(np.maximum(T_cc_te1, 1e-12))
























19:31:09 EDT Started working on Batch containing 2 tasks.                       







19:31:11 EDT Maximum FlexCredit cost: 1.051 for the whole batch.                







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












19:31:34 EDT Batch complete.                                                    



















19:31:35 EDT WARNING: Loading 'permittivity' without data; constructing a vacuum
             medium instead.                                                    







             WARNING: Loading 'permittivity' without data; constructing a vacuum
             medium instead.                                                    
























19:31:37 EDT Started working on Batch containing 2 tasks.                       







19:31:39 EDT Maximum FlexCredit cost: 1.051 for the whole batch.                







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












19:32:09 EDT Batch complete.                                                    



















19:32:10 EDT WARNING: Loading 'permittivity' without data; constructing a vacuum
             medium instead.                                                    







             WARNING: Loading 'permittivity' without data; constructing a vacuum
             medium instead.                                                    











# Reproduce paper Fig. 4c: TE0 (solid) and TE1 (dashed), with/without contour constraint
fig, ax = plt.subplots(figsize=(9, 6))

ax.plot(
    ldas_analysis * 1000, IL_std_te0, "b-", linewidth=2.0, label="Standard TO — TE0"
)
ax.plot(
    ldas_analysis * 1000,
    IL_cc_te0,
    "r-",
    linewidth=2.0,
    label="Contour-constrained — TE0",
)
ax.plot(
    ldas_analysis * 1000, IL_std_te1, "b--", linewidth=2.0, label="Standard TO — TE1"
)
ax.plot(
    ldas_analysis * 1000,
    IL_cc_te1,
    "r--",
    linewidth=2.0,
    label="Contour-constrained — TE1",
)

ax.set_xlabel("Wavelength (nm)", fontsize=12)
ax.set_ylabel("Insertion Loss (dB)", fontsize=12)
ax.invert_yaxis()
ax.set_title("O-Band Modal Demultiplexer: Insertion Loss vs. Wavelength", fontsize=13)
ax.legend(fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig("demux_il_comparison.png", dpi=150, bbox_inches="tight")
plt.show()

print("\n" + "=" * 60)
print("INSERTION LOSS SUMMARY (1260–1360 nm)")
print("=" * 60)
print(
    f"Standard TO            TE0 IL: {IL_std_te0.min():.4f} – {IL_std_te0.max():.4f} dB"
)
print(
    f"Standard TO            TE1 IL: {IL_std_te1.min():.4f} – {IL_std_te1.max():.4f} dB"
)
print(
    f"Contour-constrained    TE0 IL: {IL_cc_te0.min():.4f} – {IL_cc_te0.max():.4f} dB"
)
print(
    f"Contour-constrained    TE1 IL: {IL_cc_te1.min():.4f} – {IL_cc_te1.max():.4f} dB"
)
print()
print(
    f"TE0 improvement (paper reports ~0.010 dB): {(IL_std_te0 - IL_cc_te0).mean():+.4f} dB"
)
print(
    f"TE1 improvement (paper reports ~0.025 dB): {(IL_std_te1 - IL_cc_te1).mean():+.4f} dB"
)






















============================================================
INSERTION LOSS SUMMARY (1260–1360 nm)
============================================================
Standard TO            TE0 IL: 0.6301 – 0.7177 dB
Standard TO            TE1 IL: 0.6684 – 0.7669 dB
Contour-constrained    TE0 IL: 0.4092 – 0.5513 dB
Contour-constrained    TE1 IL: 0.6145 – 0.6876 dB

TE0 improvement (paper reports ~0.010 dB): +0.2290 dB
TE1 improvement (paper reports ~0.025 dB): +0.0528 dB










10. Device Geometry and Field Patterns


Top row: final device geometries from each branch. Bottom two rows: \(E_y\) field patterns for TE0 and TE1 input — replicates panels (a) and (b) of the paper’s Fig. 4.






fig, axes = plt.subplots(3, 2, figsize=(13, 11))

sim_std_te0_bb = get_sim(
    params_standard_seeded, beta=beta_max, input_mode=0, broadband=True, with_field=True
)
sim_cc_te0_bb = get_sim(
    params_contour_seeded, beta=beta_max, input_mode=0, broadband=True, with_field=True
)

# Geometries
sim_std_te0_bb.plot_eps(z=0.01, ax=axes[0, 0], monitor_alpha=0.0, source_alpha=0.0)
axes[0, 0].set_title("Standard TO — geometry")
sim_cc_te0_bb.plot_eps(z=0.01, ax=axes[0, 1], monitor_alpha=0.0, source_alpha=0.0)
axes[0, 1].set_title("Contour-constrained — geometry")

# TE0-in fields
batch_std["std_te0"].plot_field("field", "Ey", val="real", z=0, ax=axes[1, 0], f=freq0)
axes[1, 0].set_title("Standard TO — TE0 in (Ey, real)")
batch_cc["cc_te0"].plot_field("field", "Ey", val="real", z=0, ax=axes[1, 1], f=freq0)
axes[1, 1].set_title("Contour-constrained — TE0 in (Ey, real)")

# TE1-in fields
batch_std["std_te1"].plot_field("field", "Ey", val="real", z=0, ax=axes[2, 0], f=freq0)
axes[2, 0].set_title("Standard TO — TE1 in (Ey, real)")
batch_cc["cc_te1"].plot_field("field", "Ey", val="real", z=0, ax=axes[2, 1], f=freq0)
axes[2, 1].set_title("Contour-constrained — TE1 in (Ey, real)")

plt.tight_layout()
plt.savefig("demux_device_comparison.png", dpi=150, bbox_inches="tight")
plt.show()

























11. Conclusion


This notebook reproduced the modal demultiplexer experiment from §III-B of Wong et al. [1]. The key observations are consistent with the paper:


    

  • The contour constraint produces a visibly simpler topology (fewer small, disconnected features) than standard TO at the same iteration count.
    • 

  • Per-mode insertion-loss improvements are small but positive on average; the paper reports roughly +0.010 dB for TE0 and +0.025 dB for TE1.
    • 

  • The simplified topology is expected to be more robust to over/under-etch (process variation), although that study is not reproduced here.
    • 



Mechanistically, the contour constraint adds a gradient term that gradually erodes connected components whose integrated field intensity falls below a fraction \(T_P\) of the design’s total. For this demultiplexer we use the TE0-in forward simulation as the reference field for marking contours, since it is the dominant routing pathway.




Tuning knobs


    

  • 
T_P — threshold for marking contours (lower = more conservative, more contours preserved).
    • 

  • 
contour_weight — penalty strength (higher = faster erosion but more risk of damaging the main contours).
    • 

  • 
mask_update_freq — how often the contour mask is recomputed during Phase 2.
    • 

  • 
num_steps_phase1, num_steps_phase2, num_steps_seeded — the length of each optimization stage.
    • 







Note on iteration count


This notebook uses reduced iteration counts (num_steps_phase1=75, num_steps_phase2=25, num_steps_seeded=10) to keep the runtime short. For better final performance—particularly lower insertion loss—increase these values. The reference paper uses roughly num_steps_phase1=135 (with unconstrained_end=50, ramp_end=120), num_steps_phase2=50, and num_steps_seeded=25. The seeded-TO phase (Phase 3) benefits most from additional iterations, as it is still actively converging at 10 steps.








References


[1] J. J. Wong, R. P. Pesch, J. M. Hiesener, and S. E. Ralph, “Simplifying Photonic Topologies in Inverse Design With Contour Constraints,” IEEE Photonics Technology Letters, vol. 38, no. 2, pp. 117–120, Jan. 2026, doi: 10.1109/LPT.2025.3623046.





          

          

            Example library
          

          

          
          

          

            
            
              
              Previous Tutorial
            
            
            
            
              Next Tutorial