Waveguide Crossing with Cosine Tapers

Quentin Wach · 2026-06-25

Waveguide Crossing with Cosine Tapers

This notebook models a compact silicon waveguide crossing with cosine tapers using BEAMZ instead of Tidy3D as shown in this reference notebook. I encourage you to open up the notebooks side-by-side and compare the code and resulting plots as BEAMZ agrees with Tidy3D incredibly closely.

header_image

To achieve high integration density on a photonic chip, efficient routing of light with compact, low-loss junctions is necessary. Waveguide crossings are therefore important building blocks in high-performance photonic integrated circuits.

The convex cosine taper focuses the guided mode at the crossing center so that most power is transmitted to the through port instead of scattering into the cross ports. The design is adapted from Sujith Chandran, et al. "Beam shaping for ultra-compact waveguide crossings on monolithic silicon photonics platform," Opt. Lett. 45, 6230-6233 (2020).

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import fsolve
from shapely.geometry import Polygon as ShapelyPolygon
from shapely.geometry import box as shapely_box
from shapely.ops import unary_union

import beamz as bz
import importlib
import beamz.visual.mpl
importlib.reload(beamz.visual.mpl)

um = bz.um # unit scale in µm

plt.rcParams.update({"figure.dpi": 120})
print(f"BEAMZ version: {bz.__version__}")
BEAMZ version: 0.4.2

Simulation Setup

Define geometric parameters and materials. As in the reference design, the silicon waveguide has a thickness of 161 nm and an input width of 350 nm.

h = 0.161 * um  # waveguide thickness
w_in = 0.35 * um  # input taper width
w_out = 1.1 * um  # output taper width at the crossing center
w_m = 0.75 * um  # amplitude of the cosine function
l_t = 5.3 * um  # taper length
l_wg = 1.0 * um  # input/output straight waveguide length
n_si = 3.67
n_sio2 = 1.45

mat_si = bz.Material(permittivity=n_si**2)
mat_sio2 = bz.Material(permittivity=n_sio2**2)

The taper width is described by a cosine function $w(x)=w_m\cos(ax+b)$. We solve for $a$ and $b$ so that the taper starts at $w_\mathrm{in}/2$ and reaches $w_\mathrm{out}/2$ at the crossing junction. As in the reference notebook, one taper is generated first and the other three are made by symmetry. Before rasterization, the silicon pieces are merged into one polygon so shared boundaries do not create tiny grid seams.

# numerically solve for the cosine function that describes the taper shape
def equations(x0):
    a, b = x0
    return (
        w_m * np.cos(a * (-w_out / 2) + b) - w_out / 2,
        w_m * np.cos(a * (-w_out / 2 - l_t) + b) - w_in / 2,
    )


a, b = fsolve(equations, (0.5 / um, 2.0))

x = np.linspace(-w_out / 2 - l_t, -w_out / 2, 30)
w = w_m * np.cos(a * x + b)

# Construct one taper, then generate the remaining three by symmetry.
vertices_left = list(zip(np.r_[x, x[::-1]], np.r_[w, -w[::-1]]))
vertices_right = [(-xv, yv) for xv, yv in vertices_left]
vertices_top = [(yv, -xv) for xv, yv in vertices_left]
vertices_bottom = [(yv, xv) for xv, yv in vertices_left]
taper_vertex_sets = (vertices_left, vertices_right, vertices_top, vertices_bottom)


def _iter_polygons(geometry):
    if geometry.geom_type == "Polygon":
        return [geometry]
    if geometry.geom_type == "MultiPolygon":
        return list(geometry.geoms)
    raise ValueError(f"Expected Polygon or MultiPolygon, got {geometry.geom_type}")


def _ring_vertices(ring, z0=-h / 2):
    return [(float(xv), float(yv), z0) for xv, yv in list(ring.coords)[:-1]]


def silicon_union_2d(domain_size):
    lx, ly = domain_size
    pieces = [ShapelyPolygon(verts) for verts in taper_vertex_sets]
    pieces.extend(
        [
            shapely_box(-w_out / 2, -w_out / 2, w_out / 2, w_out / 2),
            shapely_box(-lx / 2, -w_in / 2, lx / 2, w_in / 2),
            shapely_box(-w_in / 2, -ly / 2, w_in / 2, ly / 2),
        ]
    )
    merged = unary_union(pieces)
    return merged if merged.is_valid else merged.buffer(0)


def merged_silicon_structures(domain_size, z0=-h / 2):
    return [
        bz.Polygon(
            vertices=_ring_vertices(polygon.exterior, z0),
            interiors=[_ring_vertices(ring, z0) for ring in polygon.interiors],
            material=mat_si,
            depth=h,
            z=z0,
        )
        for polygon in _iter_polygons(silicon_union_2d(domain_size))
    ]


fig, ax = plt.subplots(figsize=(5, 5))
for polygon in _iter_polygons(silicon_union_2d((2 * l_t + w_out + 2 * l_wg, 2 * l_t + w_out + 2 * l_wg))):
    x_poly, y_poly = np.asarray(polygon.exterior.coords.xy) / um
    ax.fill(x_poly, y_poly, color="tab:blue", alpha=0.65, linewidth=0)
ax.set_aspect("equal")
ax.set_xlabel("x (um)")
ax.set_ylabel("y (um)")
ax.set_title("Cosine taper crossing layout")
plt.show()
Output

BEAMZ Rasterization Note

BEAMZ rasterizes structures into a finite voxel/material grid before simulation. If overlapping or touching silicon objects are passed separately, each object can be anti-aliased independently. That can leave tiny numerical seams or asymmetries at shared boundaries, even when the analytical geometry describes the intended symmetric shape.

This notebook keeps the physical design from the reference, but adds three BEAMZ-specific safeguards before rasterization:

  • Use exact symmetric taper transforms: build the left taper once, then generate the right, top, and bottom tapers from that same reference shape.
  • Merge silicon before rasterization: union the tapers, center square, and straight waveguides into one silicon polygon so BEAMZ rasterizes one continuous region rather than several touching regions.
  • Snap the domain to the grid: make the simulation domain an integer number of grid cells so the origin and mirror axes land symmetrically on the raster grid.

These steps do not change the physical design. They make the discrete BEAMZ material grid represent the same symmetric analytical design that Tidy3D's geometry engine handles more transparently.

Set up the simulation domain, source plane, and monitors. The source and flux planes are placed in the straight waveguide sections with explicit clearance from the PML, then BEAMZ generates the actual mode source locally with ModeSolver.

lambda0 = 1.31 * um
freq0 = bz.LIGHT_SPEED / lambda0
ldas = np.linspace(1.26, 1.36, 101)
lambdas = ldas * um
freqs = bz.LIGHT_SPEED / lambdas

fwidth = freq0 / 10
run_time = 1e-12

min_steps_per_wvl = 10
grid_spec = bz.GridSpec.auto(min_steps_per_wvl=min_steps_per_wvl, wavelength=lambda0)
dx = grid_spec.resolve_resolution(max_index=n_si)


def snap_domain_to_grid(length, resolution):
    cells = int(np.ceil(length / resolution))
    # Nudge upward by one floating-point step so BEAMZ int(length / dx)
    # cannot truncate to cells - 1 because of roundoff.
    return np.nextafter(cells * resolution, np.inf), cells


Lxy_target = 2 * l_t + w_out + 2 * l_wg
Lx, nx = snap_domain_to_grid(Lxy_target, dx)
Ly, ny = snap_domain_to_grid(Lxy_target, dx)
Lz, nz = snap_domain_to_grid(1.5 * lambda0, dx)
sim_size = (Lx, Ly, Lz)

pml_layers = 12  # Tidy3D td.PML() default
pml_t = pml_layers * lambda0 / (min_steps_per_wvl * n_si)
source_x = -Lx / 2 + l_wg / 2
through_x = Lx / 2 - l_wg / 2
cross_y = Ly / 2 - l_wg / 2

# Merge the center square, four tapers, and straight waveguides into one silicon
# polygon set before rasterization. Rasterizing one merged object avoids tiny
# anti-aliased seams where separate equal-material objects overlap or touch.
structures = merged_silicon_structures((Lx, Ly), z0=-h / 2)
for structure in structures:
    structure.z = Lz / 2 - h / 2

source_time = bz.GaussianPulse(freq0=freq0, fwidth=fwidth)
mode_spec = bz.ModeSpec(num_modes=1, target_neff=3.455)
boundary_spec = bz.BoundarySpec.all_sides(
    bz.PML(
        thickness=pml_t,
        formulation="cpml",
        m=3,
        kappa_max=3.0,
        alpha_max=0.0,
    )
)

plane_z_span = 4 * h
source_plane = bz.Box(center=(source_x, 0, 0), size=(0, 4 * w_in, plane_z_span))

field_monitor = bz.FieldMonitor(
    center=(0, 0, 0),
    size=(Lx, Ly, 0),
    freqs=[freq0],
    components=("Ex", "Ey", "Ez"),
    name="field",
)

flux_monitor_through = bz.FluxMonitor(
    center=(through_x, 0, 0),
    size=(0, 4 * w_in, plane_z_span),
    freqs=freqs,
    name="flux_through",
)

flux_monitor_cross = bz.FluxMonitor(
    center=(0, cross_y, 0),
    size=(4 * w_in, 0, plane_z_span),
    freqs=freqs,
    name="flux_cross",
)

design = bz.Design(background=mat_sio2)
for structure in structures:
    design += structure

sim0 = bz.Simulation(
    domain=sim_size,
    grid_spec=grid_spec,
    design=design,
    sources=[],
    monitors=[field_monitor, flux_monitor_through, flux_monitor_cross],
    boundary_spec=boundary_spec,
    run_time=run_time,
)

print(f"dx = {sim0.resolution / um:.4f} um")
print(f"grid = {sim0.fields.permittivity.shape}, steps = {sim0.num_steps}")
print(f"target Lx/Ly = {Lxy_target / um:.4f} um; snapped Lx/Ly = {Lx / um:.4f} um ({nx} cells)")
print(f"snapped Lz = {Lz / um:.4f} um ({nz} cells)")
print(f"source/monitor x/y clearance to domain edge = {l_wg / (2 * um):.2f} um")
print(f"CPML width = {pml_t / um:.3f} um ({pml_layers} cells at the reference grid)")
print(f"source/monitor clearance to z CPML = {(Lz / 2 - pml_t - plane_z_span / 2) / um:.2f} um")

● Done: Raster cache hit (3d): 4ab95b24bed2dd13001c46cf46f58881947380e7944a4b18a911160c4e442be2.npz | load=0.19s
dx = 0.0357 um
grid = (56, 384, 384), steps = 14695
target Lx/Ly = 13.7000 um; snapped Lx/Ly = 13.7068 um (384 cells)
snapped Lz = 1.9989 um (56 cells)
source/monitor x/y clearance to domain edge = 0.50 um
CPML width = 0.428 um (12 cells at the reference grid)
source/monitor clearance to z CPML = 0.25 um

Solve the input waveguide modes on the source plane and convert the lowest-order TE-like mode into a BEAMZ mode source.

mode_solver = bz.ModeSolver(simulation=sim0, plane=source_plane, mode_spec=mode_spec, freqs=[freq0])
modes = mode_solver.solve()

fig, axes, neffs = mode_solver.plot_field_components(
    field_names=("Ey", "Ez"),
    mode_indices=(0,),
    val="abs",
    f=freq0,
    figsize=(8, 3),
    show=False,
)
display(modes.to_dataframe())
plt.show()

mode_source = mode_solver.to_source(
    mode_index=0,
    direction="+",
    source_time=source_time,
    polarization="te",
    power=1.0,
)

sim = sim0.copy(update={"sources": [mode_source]})
wavelength n eff k eff loss (dB/cm) TE (Ey) fraction wg TE fraction wg TM fraction mode area
f mode_index
2.288492e+14 0 1.31 2.343696 0.0 0.0 0.966931 0.542342 0.009508 0.099928
Output

Before running the simulation, plot the structure, source, and monitor planes to verify the setup. BEAMZ does not need a cloud job submission step; the simulation runs locally.

sim.plot(z=0, y=0, width_ratios=[1, 1.2])
plt.show()
Output

Run the local FDTD simulation. This setup now mirrors the Tidy3D reference grid density and run time, so it is significantly heavier than the earlier laptop-oriented approximation.

sim_data = sim.run(progress=True)
● JIT compiling v0.3 packed FDTD program... done!
● Progress: 100% (14695/14695 steps)

Result Visualization

After the simulation completes, first plot the field distribution at the crossing plane. A strong field focus near the junction indicates that the cosine taper is shaping the mode before it reaches the crossing center.

sim_data.plot_field(
    field_monitor_name="field",
    field_name="E",
    val="abs^2",
    f=freq0,
    vmin=0,
    vmax=3000,
    cmap="magma",
    figsize=(6, 5)
)
plt.show()
Output

Finally, quantify the crossing performance by plotting through-port insertion loss and cross-port crosstalk over the O-band wavelength range.

f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(7, 3))

T_through = sim_data["flux_through"].flux
T_cross = sim_data["flux_cross"].flux

ax1.plot(ldas, 10 * np.log10(T_through), lw=3)
ax1.set_xlabel(r"Wavelength ($\mu m$)")
ax1.set_ylabel("Transmission (dB)")
ax1.set_xlim((1.26, 1.36))
ax1.set_ylim((-0.3, -0))

ax2.plot(ldas, 10 * np.log10(T_cross / T_through), lw=3)
ax2.set_xlabel(r"Wavelength ($\mu m$)")
ax2.set_ylabel("Crosstalk (dB)")
ax2.set_xlim((1.26, 1.36))
ax2.set_ylim((-31, -25))
plt.show()
Output