BeamZWaveguide Crossing with Cosine Tapers
This notebook models a compact silicon waveguide crossing with cosine tapers using BEAMZ. It demonstrates the geometry, simulation, and analysis workflow from end to end.
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).
Related example: Tidy3D cosine-taper crossing.
import os
from pathlib import Path
import matplotlib.pyplot as plt
import numpy as np
try:
from IPython.display import display
except ImportError:
display = print
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
from beamz.simulation.backend import cuda_backend_status
test_mode = os.environ.get("BEAMZ_DOCS_TEST") == "1"
backend = "cuda_streamed"
cuda_status = cuda_backend_status()
um = bz.um # unit scale in umSimulation 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 lengthn_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 . We solve for and so that the taper starts at and reaches 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()
Keep the physical region at 13.7 × 13.7 × 1.965 µm. The source and output flux planes stay at ±6.35 µm, giving 0.5 µm clearance to the absorbing layers. Add CPML padding outside that region; use silicon's wavelength for the x/y padding and silica's wavelength for z, where the outer cells are in the cladding.
Use GridSpec.auto at 20 cells per wavelength, with grid edges at the physical-region boundaries and slab interfaces. This resolves the 161 nm silicon thickness with 10 cells. The explicit 100-million-cell budget accommodates this full 3D example.
The Gaussian center delay is 5 / (2π fwidth). Monitor data is normalized by the source spectrum. The CPML uses cubic grading, kappa_max=3, alpha_max=0, and peak conductivity 3ε₀/Δt in SI units.
lambda0 = 1.31 * um
freq0 = bz.LIGHT_SPEED / lambda0
ldas = np.linspace(1.26, 1.36, 5 if test_mode else 101)
lambdas = ldas * um
freqs = bz.LIGHT_SPEED / lambdas
fwidth = freq0 / 10
run_time = 6 / freq0 if test_mode else 1e-12
min_steps_per_wvl = 4 if test_mode else 20
pml_layers = 4 if test_mode else 12
# Keep the physical region fixed and add absorbing-layer padding around it.
Lx = Ly = 2 * l_t + w_out + 2 * l_wg
Lz = 1.5 * lambda0
physical_size = (Lx, Ly, Lz)
pml_xy = pml_layers * lambda0 / (min_steps_per_wvl * n_si)
pml_z = pml_layers * lambda0 / (min_steps_per_wvl * n_sio2)
sim_size = (Lx + 2 * pml_xy, Ly + 2 * pml_xy, Lz + 2 * pml_z)
grid_spec = bz.GridSpec.auto(
min_steps_per_wvl=min_steps_per_wvl,
wavelength=lambda0,
max_total_cells=100_000_000,
snapping_points=(
(-Lx / 2, -Ly / 2, -Lz / 2),
(Lx / 2, Ly / 2, Lz / 2),
(None, None, -h / 2),
(None, None, h / 2),
),
)
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(sim_size[:2], z0=-h / 2)
# Center the pulse at 5/(2π fwidth); offset is expressed in 1/fwidth units.
source_time = bz.GaussianPulse(
freq0=freq0,
fwidth=fwidth,
offset=0.5 if test_mode else 5.0 / (2 * np.pi),
)
mode_spec = bz.ModeSpec(num_modes=1, target_neff=3.455, polarization="te")
boundaries = tuple(
bz.PML(
edges=edges,
thickness=thickness,
formulation="cpml",
m=3,
kappa_max=3.0,
alpha_max=0.0,
)
for edges, thickness in (
(("left", "right", "bottom", "top"), pml_xy),
(("front", "back"), pml_z),
)
)
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.0, 0.0),
size=(Lx, Ly, 0.0),
freqs=[freq0],
fields=("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],
boundaries=boundaries,
run_time=run_time,
)
# Set the peak CPML conductivity in SI units after resolving the time step.
dt = float(sim0.time[1] - sim0.time[0])
boundaries = tuple(boundary.updated_copy(sigma_max=3 * bz.EPS_0 / dt) for boundary in boundaries)
sim0 = sim0.updated_copy(boundaries=boundaries)
grid_shape = sim0.grid.shape
z_edges = sim0.grid.z_edges - sim0.coordinate_offset[2]
slab_cells = np.count_nonzero((z_edges[:-1] >= -h / 2 - 1e-15) & (z_edges[1:] <= h / 2 + 1e-15))
print(f"automatic grid = {grid_shape}, cells = {np.prod(grid_shape):,}, steps = {sim0.num_steps}")
print(f"minimum spacing = {sim0.grid.minimum_spacing / bz.nm:.2f} nm; slab cells = {slab_cells}")
print(f"physical region (µm) = {np.asarray(physical_size) / um}")
print(f"domain including CPML (µm) = {np.asarray(sim_size) / um}")
print(f"CPML padding: x/y = {pml_xy / um:.3f} µm; z = {pml_z / um:.3f} µm")
print(f"source/through/cross positions (µm) = {source_x / um:.2f}, {through_x / um:.2f}, {cross_y / um:.2f}")
print(f"source/monitor clearance to x/y CPML = {l_wg / (2 * um):.2f} µm")
print(f"source/monitor clearance to z CPML = {(Lz / 2 - plane_z_span / 2) / um:.3f} µm")
assert np.isclose(source_x, -6.35 * um, rtol=1e-12, atol=1e-18)
assert np.isclose(through_x, 6.35 * um, rtol=1e-12, atol=1e-18)
assert np.isclose(cross_y, 6.35 * um, rtol=1e-12, atol=1e-18)
assert np.isclose(sim_size[0] / 2 - pml_xy - abs(source_x), l_wg / 2, rtol=1e-12, atol=1e-18)
assert np.isclose(sim_size[2] / 2 - pml_z, Lz / 2, rtol=1e-12, atol=1e-18)
assert test_mode or slab_cells >= 10
automatic grid = (792, 792, 78), cells = 48,926,592, steps = 30499 minimum spacing = 16.10 nm; slab cells = 10 physical region (µm) = [13.7 13.7 1.965] domain including CPML (µm) = [14.12833787 14.12833787 3.04913793] CPML padding: x/y = 0.214 µm; z = 0.542 µm source/through/cross positions (µm) = -6.35, 6.35, 6.35 source/monitor clearance to x/y CPML = 0.50 µm source/monitor clearance to z CPML = 0.660 µm
Solve the input waveguide modes on the source plane and convert the lowest-order TE-like mode into a BEAMZ mode source.
mode_source_request = bz.ModeSource(
center=source_plane.center,
size=source_plane.size,
direction="+",
source_time=source_time,
mode_spec=mode_spec,
power=1.0,
)
modes = mode_source_request.solve_modes(sim0, freqs=[freq0])
from beamz.analysis import mode_data_to_dataframe
from beamz.analysis.plotting import plot_mode_field_components
fig, axes, neffs = plot_mode_field_components(modes,
field_names=("Ey", "Ez"),
mode_indices=(0,),
val="abs",
f=freq0,
figsize=(8, 3),
show=False,
)
display(mode_data_to_dataframe(modes))
plt.show()
# Attach the same immutable ModeSource config after inspection. Compilation will
# lower it into a native launch plan without mutating this object.
mode_source = mode_source_request
sim = sim0.updated_copy(sources=(mode_source,))| f | mode_index | wavelength | n eff | k eff | loss (dB/cm) | mode area |
|---|---|---|---|---|---|---|
| 2.288492e+14 | 0 | 1.31 | 2.32836 | 0.0 | 0.0 | 0.096229 |
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()Run the crossing with the automatic mesh at 20 cells per wavelength for the full 1 ps. ModeSource calibrates its emitted power as part of launch compilation, and the detached flux results apply the source waveform normalization automatically. No reference or input monitor is required.
sim_data = sim.run(progress=True, backend=backend)* Compiling simulation... Simulation runtime: 275.06 s GCUPS: 5.425
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()Finally, quantify the crossing performance over the O-band wavelength range. FluxMonitor.flux already includes the source waveform normalization and is reported in watts. The transmission denominator is SimulationResults.launched_power(): the solved source mode integrated with the mode solver's exact component-staggered transverse metric. The mode-source compiler uses that value to calibrate the launch to the requested 1 W, without adding an input monitor or reference simulation. Crosstalk is the cross-port power relative to through-port power. We plot 10 log10(T), which is negative for a passive device; positive insertion loss would instead use -10 log10(T).
f, (ax1, ax2) = plt.subplots(1, 2, tight_layout=True, figsize=(7, 3))
flux_through = sim_data["flux_through"].flux
flux_cross = sim_data["flux_cross"].flux
source_power = sim_data.launched_power(source=0)
def safe_power_ratio(numerator, denominator, floor=1e-30):
numerator = np.asarray(numerator, dtype=float)
denominator = np.asarray(denominator, dtype=float)
ratio = np.zeros(np.broadcast_shapes(numerator.shape, denominator.shape), dtype=float)
return np.divide(numerator, denominator, out=ratio, where=np.abs(denominator) >= floor)
T_through = safe_power_ratio(flux_through, source_power)
T_cross = safe_power_ratio(flux_cross, source_power)
transmission_db = 10 * np.log10(np.maximum(T_through, 1e-30))
crosstalk_db = 10 * np.log10(np.maximum(safe_power_ratio(T_cross, T_through), 1e-30))
ax1.plot(ldas, transmission_db, 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, crosstalk_db, 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()