Waveguide Crossing Yield Analysis

Dense photonic circuits route hundreds of waveguide crossings, so a crossing must not only perform well at its nominal design, it must keep performing across the spread of silicon thickness and linewidth that a real fabrication process delivers. This example predicts the fabrication yield of an adiabatic waveguide crossing on the SiEPIC OpenEBL platform using the response-surface methodology of experimental design [1]:

1. Simulate the crossing with full 3D FDTD at a small set of process corners and samples.

2. Fit smooth polynomial response surfaces to the results, and check them against held-out simulations the fit has never seen.

3. Sample hundreds of virtual chips from a hierarchical wafer and lot variation model [2], evaluating each chip through the fitted surfaces in microseconds.

The full-wave cost is fixed by step 1 (eleven simulations here), independent of how many chips step 3 samples. This is what makes statistical tolerance analysis affordable: the simulator explores the process space once, and the statistics are then free.

Reference

1. Box, George E. P., and Kenneth B. Wilson. “On the experimental attainment of optimum conditions.” Journal of the Royal Statistical Society: Series B 1951 13 (1), 1-45, doi: 10.1111/j.2517-6161.1951.tb00067.x.

2. Lu, Zeqin, et al. “Performance prediction for silicon photonics integrated circuits with layout-dependent correlated manufacturing variability.” Optics Express 2017 25 (9), 9712-9733, doi: 10.1364/OE.25.009712.

3. Ma, Yangjin, et al. “Ultralow loss single layer submicron silicon waveguide crossing for SOI optical interconnect.” Optics Express 2013 21 (24), 29374-29382, doi: 10.1364/OE.21.029374.

Technology and variation axes

The two dominant process variations on a silicon photonics platform are the device-layer thickness (a wafer-scale property of the SOI stock) and the waveguide linewidth, or critical dimension (CD), set by lithography and etch. Both are exposed as parameters of the SiEPIC OpenEBL technology: si_thickness sets the extruded silicon thickness, and si_mask_dilation moves every silicon mask edge outward by the given amount, so a dilation of \(d\) changes every width by \(2d\), exactly like a litho or etch bias.

One setup detail matters here: the crossing uses 430 nm access waveguides, so a port specification for that width is created (by copying and adjusting the stock 500 nm TE port) and passed directly to the component.

import numpy as np
import photonforge as pf
import siepic_forge as siepic
from scipy.interpolate import make_interp_spline
from matplotlib import pyplot as plt
from photonforge.live_viewer import LiveViewer

viewer = LiveViewer()

# SiEPIC OpenEBL technology at its nominal stack (220 nm silicon, no bias)
tech = siepic.ebeam()

# Port specification for the 430 nm access waveguides, made by copying and
# adjusting the stock 500 nm TE port. It is passed directly to the component,
# so it does not need to be registered in the technology.
wg_width = 0.43
port_spec = tech.ports["TE_1550_500"].copy()
port_spec.path_profiles = [(wg_width, 0.0, (1, 0))]
port_spec.description = "Strip TE 1550 nm, w=430 nm"

pf.config.default_technology = tech

# Analysis band: 11 wavelengths across the C band
wavelengths = np.linspace(1.53, 1.57, 11)
freqs = pf.C_0 / wavelengths

11:22:07 Eastern Daylight Time WARNING: The material-library variant            
                               'Palik_Lossless' is deprecated and maps to       
                               'Palik_LowLoss' because it contains a tiny fitted
                               loss despite its name. Use 'Palik_NoLoss' where  
                               available for a zero-loss Palik model.           

LiveViewer started at http://localhost:52751

Adiabatic crossing component

The crossing consists of four identical arms joined at the center. Each arm expands the 430 nm access waveguide through a cubic-spline width profile up to a wide multimode section at the junction, where the mode is strongly confined and expands slowly enough to cross the perpendicular arms with little scattering, the approach that established ultralow crossing loss on submicron SOI [3]. The spline knots (default_widths) were chosen so the expansion stays adiabatic over the C band.

The Tidy3DModel is attached inside the parametric function. The declared port_symmetries (both mirror planes and inversion) tell the model that a single FDTD excitation determines the response of symmetry-related inputs, which cuts the number of simulations per geometry.

# Spline knot widths from the arm tip (port side) to the junction. The wide
# 2.15 um section sits next to the center; the final 0.5 um value is the
# width at the exact center, where the perpendicular arms overlap.
default_widths = (
    0.5, 0.6, 0.95, 1.32, 1.44, 1.46, 1.466,
    1.52, 1.58, 1.62, 1.76, 2.15, 0.5,
)


@pf.parametric_component(name_prefix="PLUS_CROSSING")
def adiabatic_plus_crossing(
    *,
    port_spec,
    arm_length=4.5,
    widths=default_widths,
    lead_length=1.0,
):
    if isinstance(port_spec, str):
        port_spec = pf.config.default_technology.ports[port_spec]

    wg_width, _ = port_spec.path_profile_for("Si")

    # Number of points used to discretize the cubic spline profile.
    num_points = int(arm_length / pf.config.tolerance)

    # Snap the port coordinate to the grid to avoid gaps/overlaps.
    xp = pf.snap_to_grid(arm_length + lead_length)

    # Cubic spline through the width knots (reversed so that 0 = center).
    coords = np.linspace(0, arm_length, len(widths))
    spline = make_interp_spline(coords, widths[::-1], k=3)

    # Pre-compute widths and positions from the interpolation.
    coords = np.linspace(arm_length, 0, num_points)
    w = spline(coords)

    # East arm: lead-in taper from the port, then the spline profile.
    arm_e = pf.Path((xp, 0), wg_width)
    arm_e.segment((coords[0], 0), width=(w[0], "smooth"))
    for x, wi in zip(coords[1:], w[1:]):
        arm_e.segment((x, 0), wi)

    # The other three arms are rotated copies of the east arm.
    arm_w = arm_e.copy().rotate(180)
    arm_n = arm_e.copy().rotate(90)
    arm_s = arm_e.copy().rotate(270)

    # Merge the four arms into a single silicon region.
    c = pf.Component()
    c.add("Si", *pf.boolean([arm_e, arm_w], [arm_n, arm_s], "+"))

    c.add_port(c.detect_ports([port_spec]))
    assert len(c.ports) == 4, "Port detection failed: expected exactly 4 ports."

    # Detected port order: P0 = West, P1 = South, P2 = North, P3 = East.
    c.add_model(
        pf.Tidy3DModel(
            grid_spec=12,
            port_symmetries=[
                ("P0", "P2", "P1", "P3"),  # symmetry about x-axis
                ("P3", "P1", "P2", "P0"),  # symmetry about y-axis
                ("P3", "P2", "P1", "P0"),  # inversion symmetry
            ],
        ),
        "Tidy3DModel",
    )
    return c


crossing = adiabatic_plus_crossing(port_spec=port_spec)
viewer(crossing)

Nominal performance

A crossing is judged on three quantities: the through transmission (P0 to P3), the crosstalk scattered into the perpendicular arms (P0 to P1 and P2), and the back-reflection. We compute the full scattering matrix over the band and adopt the following acceptance criteria, evaluated at the worst wavelength: transmission above 98% and crosstalk below -30 dB.

# Full S-matrix of the nominal crossing (one FDTD task after symmetries)
s_nominal = crossing.s_matrix(freqs)


def crossing_metrics(s):
    """Through power, worst crosstalk (dB) and reflection (dB) vs wavelength."""
    # Through path: west input to east output
    thru = np.abs(s[("P0@0", "P3@0")]) ** 2
    # Crosstalk: worst of the two perpendicular outputs
    xt = np.maximum(
        np.abs(s[("P0@0", "P1@0")]) ** 2,
        np.abs(s[("P0@0", "P2@0")]) ** 2,
    )
    # Back-reflection into the input port
    refl = np.abs(s[("P0@0", "P0@0")]) ** 2
    return thru, 10 * np.log10(xt), 10 * np.log10(refl)


t_nom, xt_nom, r_nom = crossing_metrics(s_nominal)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 3.6), tight_layout=True)
ax1.plot(wavelengths, 100 * t_nom, "o-")
ax1.axhline(98, color="tab:red", ls="--", label="98% criterion")
ax1.set_xlabel("wavelength (um)")
ax1.set_ylabel("transmission (%)")
ax1.legend()
ax1.grid(alpha=0.3)
ax2.plot(wavelengths, xt_nom, "o-", label="crosstalk")
ax2.plot(wavelengths, r_nom, "s--", label="reflection")
ax2.axhline(-30, color="tab:red", ls="--", label="-30 dB criterion")
ax2.set_xlabel("wavelength (um)")
ax2.set_ylabel("power (dB)")
ax2.legend()
ax2.grid(alpha=0.3)
plt.show()

print(f"worst-band transmission: {100 * t_nom.min():.2f}%")
print(f"worst-band crosstalk: {xt_nom.max():.1f} dB")

11:22:13 Eastern Daylight Time Loading simulation from local cache. View cached
                               task using web UI at
                               'https://tidy3d.simulation.cloud/workbench?taskId
                               =fdve-0e4329f8-263c-4425-a949-679caa9ffa09'.

                               Loading simulation from local cache. View cached
                               task using web UI at
                               'https://tidy3d.simulation.cloud/workbench?taskId
                               =fdve-a25098c6-24a5-4ab7-b001-2a75fdf95c38'.

worst-band transmission: 98.48%
worst-band crosstalk: -58.1 dB

The nominal crossing transmits above 98.4% everywhere in the band, with the minimum at the red edge where the adiabatic expansion is slowly running out of margin. Crosstalk and reflection sit below -55 dB and -40 dB, far from their criterion. The 0.47 point transmission margin at the band edge is the quantity the process variation will eat into, so the yield question is a real one.

Process corner simulations

We now map how the metrics move across the process space. The monte_carlo.s_matrix function automates the whole batch: the two process deviations are declared as random variables of the technology, using the platform budgets as normal distributions (\(\sigma_t = 4\) nm for thickness and \(\sigma_w = 3.3\) nm for linewidth), and a single call evaluates

• the four corner cases of the process box, whose edges are placed at 2.5 standard deviations by corner_coverage, and

• six random samples inside the box (random_samples).

Technology updates, component regeneration, parallel upload, and collection are all handled by the function. The four box corners alone cannot separate the two curvature terms of a quadratic surface, so the first four random samples complete the training set, and the last two are set aside as hold-outs that never enter the fit. One FDTD input suffices per case because the four-fold symmetry of the crossing makes all inputs equivalent (the nominal run, which computed two independent inputs, confirms this equivalence to numerical precision).

thickness_var = pf.monte_carlo.RandomVariable(
    "si_thickness", tech, value=0.220, stdev=0.004
)
# The mask dilation moves each silicon edge by half the linewidth change,
# so a 3.3 nm linewidth sigma is a 1.65 nm dilation sigma (values in um).
dilation_var = pf.monte_carlo.RandomVariable(
    "si_mask_dilation", tech, value=0.0, stdev=0.00165
)

_, mc_results = pf.monte_carlo.s_matrix(
    crossing,
    freqs,
    thickness_var,
    dilation_var,
    corner_samples=-1,  # all 4 corners of the process box
    random_samples=6,  # random samples inside the box
    corner_coverage=0.9876,  # box edges at 2.5 standard deviations
    random_seed=0,
    model_kwargs={"inputs": ["P0"]},
)

# Each result is (thickness, dilation, s_matrix). Convert the sampled
# absolute values to deviations in nm for bookkeeping.
t_corner, xt_corner = {}, {}
for thickness, dilation, s in mc_results:
    key = (round(1e3 * (thickness - 0.220), 1), round(2e3 * dilation, 1))
    thru, xt, _ = crossing_metrics(s)
    t_corner[key], xt_corner[key] = thru, xt

# The 4 corners and the first 4 samples train the fit; the last 2 samples
# are reserved as hold-outs.
doe = list(t_corner)
doe_train, doe_holdout = doe[:8], doe[8:]

print(" dt (nm)  dw (nm)   worst T (%)   worst XT (dB)")
for dt, dw in doe:
    print(f"  {dt:+5.1f}   {dw:+5.1f}     {100 * t_corner[(dt, dw)].min():6.2f}       "
          f"{xt_corner[(dt, dw)].max():6.1f}")

Uploading task 'P0@0…'
Uploading task 'P0@0…'
Uploading task 'P0@0…'
Uploading task 'P0@0…'
Uploading task 'P0@0…'
Uploading task 'P0@0…'
Uploading task 'P0@0…'
Uploading task 'P0@0…'
Uploading task 'P0@0…'
Uploading task 'P0@0…'
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-88d7a0f8-0d2c-4085-bc31-430187a41b25
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-bd217112-e1a2-4b3c-a39a-5641a5492c04
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-a29b52e7-6026-4a94-8158-9cf908496e1a
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-024ae7ae-8c8d-472f-868a-29bc0615844c
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-2394966c-36ed-4615-ac6c-14b86963acf6
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-4a9dceb9-03ab-4300-892b-f7fce87d4863
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-69f453c1-4f3e-40ec-b2ed-55278f29868f
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-bfc80ae1-1652-4fc6-96e3-a514ec4eab8e
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-4fb58a39-9171-405a-95bc-3b1121b5be4c
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-7d5eb6d4-e389-4cca-b3d9-7681c79f0253
Downloading data from 'P0@0'…
Downloading data from 'P0@0'…
Downloading data from 'P0@0'…
Downloading data from 'P0@0'…
Downloading data from 'P0@0'…
Downloading data from 'P0@0'…
Downloading data from 'P0@0'…
Downloading data from 'P0@0'…
Downloading data from 'P0@0'…
Downloading data from 'P0@0'…
 dt (nm)  dw (nm)   worst T (%)   worst XT (dB)
  -10.0    -8.3      97.35        -53.0
  +10.0    -8.3      98.34        -58.8
  -10.0    +8.3      98.68        -57.8
  +10.0    +8.3      99.34        -53.1
   +1.2    +6.5      98.94        -58.2
   -1.0    +0.5      98.49        -57.9
   -1.9    -5.6      98.02        -54.7
   +7.3    -0.7      98.76        -60.1
   +3.7    -2.5      98.47        -59.4
   -4.9    +2.0      98.49        -56.9

The worst-band transmission is lowest at the thin and narrow corner and highest at the thick and wide one: more silicon always helps this crossing, and the 98% criterion boundary lies inside the process box. Crosstalk stays below -50 dB at every sampled point. Since that is more than 20 dB inside the -30 dB criterion, crosstalk cannot drive yield and we treat it as bounded by its worst sampled value instead of fitting a surface to it (at these levels the FDTD result is dominated by its numerical noise floor, so a fitted surface would model noise, not physics).

Transmission response surface

The transmission at each wavelength is fitted as a quadratic in the two deviations, anchored exactly at the nominal:

\[T(\lambda;\, \Delta t, \Delta w) = T_0(\lambda) + c_1 \Delta t + c_2 \Delta w + c_3 \Delta t^2 + c_4 \Delta w^2 + c_5 \Delta t \, \Delta w\]

The eight training cases determine the five coefficients per wavelength by least squares. The purpose of the plot below is to test the accuracy of the fitted surface: the two hold-out samples were never used in the fit, so comparing their FDTD spectra (markers) with the surface’s predictions at the same deviations (lines) measures the model’s true prediction error.

def q_basis(dt, dw):
    """Quadratic basis terms (dt, dw, dt^2, dw^2, dt*dw) for the fit."""
    dt, dw = np.asarray(dt, float), np.asarray(dw, float)
    return np.stack([dt, dw, dt**2, dw**2, dt * dw], axis=-1)


# Fit the deviation from the nominal spectrum: 8 training cases determine
# the 5 coefficients per wavelength by least squares.
a_train = q_basis([p[0] for p in doe_train], [p[1] for p in doe_train])
y_train = np.stack([t_corner[p] for p in doe_train]) - t_nom
t_coef = np.linalg.lstsq(a_train, y_train, rcond=None)[0]


def t_model(dt, dw):
    """Predicted transmission spectrum at deviation (dt, dw) in nm."""
    return t_nom + q_basis(dt, dw) @ t_coef


# The two hold-outs measure the model's true prediction error.
holdout_err = max(
    np.max(np.abs(t_model(dt, dw) - t_corner[(dt, dw)])) for dt, dw in doe_holdout
)

fig, ax = plt.subplots(figsize=(7.5, 4.2), tight_layout=True)
for (dt, dw), color in zip(doe_holdout, ("tab:blue", "tab:orange")):
    # Response-surface prediction at the held-out deviation
    ax.plot(wavelengths, 100 * t_model(dt, dw), color=color, lw=2.5,
            label=f"fit at ({dt:+.1f}, {dw:+.1f}) nm")
    # FDTD simulation at the same deviation, never used in the fit
    ax.plot(wavelengths, 100 * t_corner[(dt, dw)], "o", color=color,
            ms=7, ls="none")
ax.plot([], [], "ko", ms=7, label="FDTD (not used in the fit)")
ax.set_xlabel("wavelength (um)")
ax.set_ylabel("transmission (%)")
ax.legend()
ax.grid(alpha=0.3)
plt.show()

print(f"largest hold-out error: {100 * holdout_err:.2f} transmission points")

largest hold-out error: 0.09 transmission points

The held-out simulations land on the predicted curves to a fraction of a transmission point, an order of magnitude tighter than the margins that decide pass or fail, so the surface can stand in for the simulator across the whole box.

Wafer and lot variation analysis

Real variation is hierarchical: silicon thickness and CD differ from lot to lot, from wafer to wafer within a lot, and across each wafer [2]. We use a compact model of that hierarchy with the same total budgets as the corner variables (\(\sigma_t = 4\) nm and \(\sigma_w = 3.3\) nm), split between a lot offset, a wafer offset, a systematic across-wafer profile (radial bow for thickness, a linear tilt for CD), and a per-chip random term. Every chip is then a \((\Delta t, \Delta w)\) pair evaluated through the fitted surface, and its worst-band transmission decides pass or fail.

rng = np.random.default_rng(7)
n_lots, n_wafers, n_chips = 3, 5, 40
wafer_radius = 100.0  # mm

# Fixed chip sample plan (same die positions measured on every wafer)
angle = rng.uniform(0, 2 * np.pi, n_chips)
radius = wafer_radius * 0.85 * np.sqrt(rng.uniform(0, 1, n_chips))
site_x, site_y = radius * np.cos(angle), radius * np.sin(angle)
r2 = (site_x**2 + site_y**2) / wafer_radius**2  # normalized radius squared

records = []
for lot in range(n_lots):
    # Lot offsets: every wafer in the lot shares them
    lot_t, lot_w = rng.normal(0, 2.4), rng.normal(0, 1.7)
    for wafer in range(n_wafers):
        # Wafer offsets on top of the lot offsets
        wafer_t = lot_t + rng.normal(0, 2.4)
        wafer_w = lot_w + rng.normal(0, 1.7)
        # Chip values: systematic across-wafer profile + per-chip random term
        chip_t = wafer_t - 3.0 * (r2 - 0.5) + rng.normal(0, 1.5, n_chips)
        chip_w = wafer_w + 1.3 * site_x / wafer_radius + rng.normal(0, 2.0, n_chips)
        if lot == 0 and wafer == 0:
            # Keep the first wafer for the wafer-map figure below
            example_wafer = (wafer_t, wafer_w, chip_t.copy(), chip_w.copy())
        for dt, dw in zip(chip_t, chip_w):
            # Evaluate inside the trained box (rare outliers are clamped)
            dtc, dwc = np.clip(dt, -10.0, 10.0), np.clip(dw, -8.25, 8.25)
            records.append((lot, wafer, dt, dw, t_model(dtc, dwc).min()))

lot_idx = np.array([r[0] for r in records])
wafer_idx = np.array([r[1] for r in records])
chip_dt = np.array([r[2] for r in records])
chip_dw = np.array([r[3] for r in records])
t_worst = np.array([r[4] for r in records])

n = len(records)
passed = t_worst > 0.98
yield_est = passed.mean()

# Wilson 95% confidence interval on the yield estimate
z = 1.96
center = (yield_est + z**2 / (2 * n)) / (1 + z**2 / n)
half = z * np.sqrt(yield_est * (1 - yield_est) / n + z**2 / (4 * n**2)) / (1 + z**2 / n)

print(f"chips: {n}")
print(f"yield (worst-band T > 98%): {100 * yield_est:.1f}%  "
      f"(95% CI {100 * (center - half):.1f}-{100 * (center + half):.1f}%)")

chips: 600
yield (worst-band T > 98%): 91.8%  (95% CI 89.4-93.8%)

The hierarchy is easiest to see on a single wafer. The maps below show the thickness and linewidth deviations across the first sampled wafer: the smooth background combines the lot and wafer offsets with the systematic across-wafer profiles (a radial bow for thickness and a linear tilt for linewidth), and the dots are the sampled chip sites, which add the per-chip random term on top of the smooth map.

wafer_t0, wafer_w0, chips_t0, chips_w0 = example_wafer

# Systematic wafer-scale profiles evaluated on a grid covering the wafer
axis_mm = np.linspace(-wafer_radius, wafer_radius, 201)
gx, gy = np.meshgrid(axis_mm, axis_mm)
inside = gx**2 + gy**2 <= wafer_radius**2
rr2 = (gx**2 + gy**2) / wafer_radius**2
map_t = np.where(inside, wafer_t0 - 3.0 * (rr2 - 0.5), np.nan)
map_w = np.where(inside, wafer_w0 + 1.3 * gx / wafer_radius, np.nan)

fig, axes = plt.subplots(1, 2, figsize=(10.5, 4.4), tight_layout=True)
for ax, field, chips, label in (
    (axes[0], map_t, chips_t0, "thickness deviation (nm)"),
    (axes[1], map_w, chips_w0, "linewidth deviation (nm)"),
):
    # Shared color scale between the smooth map and the chip sites
    vmin = min(np.nanmin(field), chips.min())
    vmax = max(np.nanmax(field), chips.max())
    mesh = ax.pcolormesh(gx, gy, field, cmap="RdBu_r", vmin=vmin, vmax=vmax,
                         shading="auto")
    ax.scatter(site_x, site_y, c=chips, cmap="RdBu_r", vmin=vmin, vmax=vmax,
               s=36, edgecolors="black", linewidths=0.5)
    plt.colorbar(mesh, ax=ax, label=label)
    ax.set_aspect("equal")
    ax.set_xlabel("x (mm)")
    ax.set_ylabel("y (mm)")
axes[0].set_title("silicon thickness across one wafer")
axes[1].set_title("linewidth across one wafer")
plt.show()

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10.5, 4), tight_layout=True)

# One box per wafer, colored by lot: exposes which hierarchy level moves
data, positions, box_colors = [], [], []
cmap = plt.get_cmap("tab10")
for lot in range(n_lots):
    for wafer in range(n_wafers):
        sel = (lot_idx == lot) & (wafer_idx == wafer)
        data.append(100 * t_worst[sel])
        positions.append(lot * (n_wafers + 1.5) + wafer)
        box_colors.append(cmap(lot))
boxes = ax1.boxplot(data, positions=positions, widths=0.7, patch_artist=True)
for patch, color in zip(boxes["boxes"], box_colors):
    patch.set_facecolor(color)
ax1.axhline(98, color="tab:red", ls="--")
ax1.set_xticks([lot * (n_wafers + 1.5) + 2 for lot in range(n_lots)])
ax1.set_xticklabels([f"lot {lot + 1}" for lot in range(n_lots)])
ax1.set_ylabel("worst-band transmission (%)")
ax1.set_title("per-wafer distributions, grouped by lot")

# Population histogram against the acceptance criterion
ax2.hist(100 * t_worst, bins=40)
ax2.axvline(98, color="tab:red", ls="--", label="98% criterion")
ax2.set_xlabel("worst-band transmission (%)")
ax2.set_ylabel("chips")
ax2.set_title(f"{n} chips, yield {100 * yield_est:.1f}%")
ax2.legend()
plt.show()

The box plots expose the failure structure: whole lots move together, because the lot-level thickness and CD offsets are shared by every chip in the lot. A lot that draws thin silicon together with narrow CD straddles the criterion and contributes most of the failures, while favorable lots yield essentially 100%. For a fab this is actionable: screening incoming wafers on device-layer thickness catches most failing material before any processing.

The same data can be viewed directly in the process plane. The map below shows the worst-band transmission predicted by the response surface over the whole deviation box, the 98% criterion boundary, and every sampled chip.

# Evaluate the response surface on a grid covering the process box
dt_grid = np.linspace(-10, 10, 121)
dw_grid = np.linspace(-8.25, 8.25, 121)
dt_mesh, dw_mesh = np.meshgrid(dt_grid, dw_grid, indexing="ij")
# Worst-band value at every grid point: minimum of the spectrum
t_map = t_model(dt_mesh.ravel()[:, None], dw_mesh.ravel()[:, None])
t_map = t_map.reshape(dt_mesh.shape + (-1,)).min(axis=-1)

fig, ax = plt.subplots(figsize=(6.8, 5), tight_layout=True)
contours = ax.contourf(dt_mesh, dw_mesh, 100 * t_map, 21, cmap="RdYlGn")
plt.colorbar(contours, ax=ax, label="worst-band transmission (%)")
# The acceptance boundary in the process plane
boundary = ax.contour(dt_mesh, dw_mesh, 100 * t_map, levels=[98.0],
                      colors="k", linewidths=2)
ax.clabel(boundary, fmt="98%%")
# Every sampled chip at its process deviation
ax.plot(chip_dt, chip_dw, "k.", ms=2.5, alpha=0.35)
ax.set_xlim(-10, 10)
ax.set_ylim(-8.25, 8.25)
ax.set_xlabel("silicon thickness deviation (nm)")
ax.set_ylabel("linewidth deviation (nm)")
ax.set_title("Sampled chips over the process window")
plt.show()

Reading the surface as a design tool

The fitted coefficients also quantify how to buy margin. At the worst wavelength the surface slopes are about +0.05 transmission points per nm of silicon thickness and +0.07 points per nm of linewidth: this crossing always improves with more silicon. Thickness is not a designer’s knob, but linewidth is, since the drawn layout can be biased. Widening the drawn profile by 5 to 10 nm would roughly double the worst-band margin at zero cost, moving the criterion boundary mostly out of the process box. Re-centering the transmission peak (currently near 1.538 um) onto the band center by retuning the spline profile would have a similar effect. Both changes can be evaluated through the same corner-plus-surface workflow before committing a new layout.

The method shown here generalizes to any component whose response varies smoothly across the process box: simulate a designed set of corners once, prove the surface on held-out corners, then sample the process for free. Resonant devices need one adaptation, fitting physical parameters such as coupling and effective index instead of raw spectra, but the corner budget and the statistical machinery stay the same.