BeamZ1×4 MMI power splitter in BeamZ
This BeamZ-native workflow recreates the geometry, source, monitors, and spectral analysis of Tidy3D's 1×4 MMI power-splitter example. The device is a 3D Si–GaN–Si slot-waveguide MMI at 1550 nm.
The normal notebook uses an automatic rectilinear mesh with a 16-points-per-wavelength target and acquires 401 wavelengths from 1500 to 1600 nm (0.25 nm spacing) in every flux and mode monitor. These are directly sampled DFT values, rather than interpolated curves. The source and mode monitors select TM-like modes, with five source frequency profiles for broadband excitation. Denser wavelength sampling does not replace spatial or time-domain convergence checks.
Field-decay and monitor-stability checks can stop the run early. A 20-minute wall-clock budget stops execution at a chunk boundary and retains acquired measurements; final result construction can take the total time beyond that budget. Spectra are labeled provisional if the run stops before convergence. BEAMZ_DOCS_TEST=1 activates a small three-wavelength configuration for automated validation.
Saved execution: wall_time_limit after 224,000 steps (11.032 ps). The saved spectra are provisional. Final field decay was 1.86e-12 (threshold 1e-5), monitor change was 7.29e-08 (threshold 1e-3), and remaining per-term source activity was 0.00754 (source-off threshold 1e-6). The 1200-second execution budget retained the acquired measurements.
import os
import matplotlib.pyplot as plt
import numpy as np
from IPython.display import display
import beamz as bz
from beamz.analysis import mode_data_to_dataframe
from beamz.analysis.plotting import plot_mode_field_components
test_mode = os.environ.get("BEAMZ_DOCS_TEST") == "1"
execution_backend = "jax" if test_mode else "cuda_streamed"
um = bz.um1. Materials and design parameters
The parameters match the reference example: a 5 µm-wide MMI, 11.2 µm initial length, 0.4 µm access waveguides, and a 0.1 µm GaN slot between 0.3 µm silicon layers.
n_si = 3.48
n_gan = 2.305
n_sio2 = 1.444
mat_si = bz.Material(permittivity=n_si**2)
mat_gan = bz.Material(permittivity=n_gan**2)
mat_sio2 = bz.Material(permittivity=n_sio2**2)
W_MMI = 5.0 * um
L_MMI = 11.2 * um
g1 = 0.9 * um
W1 = 0.4 * um
W2 = 0.8 * um
L1 = 2.0 * um
L2 = 5.0 * um
H_SI = 0.3 * um
H_GAN = 0.1 * um
lambda0 = 1.55 * um
freq0 = bz.LIGHT_SPEED / lambda0
wavelength_samples = 3 if test_mode else 401
wavelengths = np.linspace(1.50, 1.60, wavelength_samples) * um
freqs = bz.LIGHT_SPEED / wavelengths
fwidth = 0.5 * (np.max(freqs) - np.min(freqs))
buffer_x = 1.0 * um
buffer_y = 1.5 * um
slot_stack_thickness = 2 * H_SI + H_GAN
domain_depth = 5 * (2 * H_GAN + H_SI)
output_centers = np.array([-1.5, -0.5, 0.5, 1.5]) * (g1 + W1)
print("Output centerlines (µm):", output_centers / um)Output centerlines (µm): [-1.95 -0.65 0.65 1.95]
2. Build a parameterized 3D MMI simulation
make_sim() is deliberately analogous to the reference notebook. It returns an immutable BeamZ Simulation and the named monitor geometry used later for analysis.
API opportunity: a public native MMI builder with named optical ports could replace the explicit rectangles and taper polygons below while preserving this transparent workflow.
def trapezoid(x0, x1, width0, width1, center_y=0.0):
"""Return a longitudinal taper polygon in the xy plane."""
return (
(x0, center_y - 0.5 * width0),
(x1, center_y - 0.5 * width1),
(x1, center_y + 0.5 * width1),
(x0, center_y + 0.5 * width0),
)
def make_sim(L_mmi, W_mmi, *, include_field=True):
domain_length = L1 + L_mmi + L2 + 2 * buffer_x
domain_width = W_mmi + 2 * buffer_y
domain = (domain_length, domain_width, domain_depth)
x_left = -0.5 * domain_length
x_taper_in = x_left + buffer_x
x_mmi_in = x_taper_in + L1
x_mmi_out = x_mmi_in + L_mmi
x_taper_out = x_mmi_out + L2
# Si is continuous above and below the GaN slot. The GaN polygons are
# added second, so the slot overrides the silicon material in its slab.
silicon_polygons = [
((x_left, -W1 / 2), (x_taper_in, -W1 / 2),
(x_taper_in, W1 / 2), (x_left, W1 / 2)),
trapezoid(x_taper_in, x_mmi_in, W1, W2),
((x_mmi_in, -W_mmi / 2), (x_mmi_out, -W_mmi / 2),
(x_mmi_out, W_mmi / 2), (x_mmi_in, W_mmi / 2)),
]
for center_y in output_centers:
silicon_polygons.extend(
[
trapezoid(x_mmi_out, x_taper_out, W2, W1, center_y),
((x_taper_out, center_y - W1 / 2), (0.5 * domain_length, center_y - W1 / 2),
(0.5 * domain_length, center_y + W1 / 2), (x_taper_out, center_y + W1 / 2)),
]
)
design = bz.Design(background=mat_sio2)
for vertices in silicon_polygons:
design += bz.Polygon(
vertices=vertices, material=mat_si, depth=slot_stack_thickness,
z=-0.5 * slot_stack_thickness, color="#d81b60",
)
for vertices in silicon_polygons:
design += bz.Polygon(
vertices=vertices, material=mat_gan, depth=H_GAN, z=-0.5 * H_GAN,
color="#f6c85f",
)
steps_per_wavelength = 4 if test_mode else 16
grid_spec = bz.GridSpec.auto(
wavelength=lambda0, min_steps_per_wvl=steps_per_wavelength, max_scale=1.3,
)
# Keep the physical PML thickness and port apertures fixed as the mesh changes.
pml_thickness = lambda0 / n_si
source_time = bz.GaussianPulse(
freq0=freq0, fwidth=fwidth, offset=0.5 if test_mode else 4.0,
)
mode_spec = bz.ModeSpec(
num_modes=2, mode_index=0, target_neff=n_si, polarization="tm", num_freqs=5,
)
source_x = x_left + 0.5 * buffer_x
output_x = 0.5 * domain_length - 0.5 * buffer_x
plane_size = (0.0, 5 * W1, domain_depth - 2 * pml_thickness)
source = bz.ModeSource(
center=(source_x, 0.0, 0.0), size=plane_size, direction="+",
source_time=source_time, mode_spec=mode_spec,
)
flux_monitors = tuple(
bz.FluxMonitor(
center=(output_x, center_y, 0.0), size=(0.0, 3 * W1, domain_depth - 2 * pml_thickness),
freqs=freqs, name=f"flux_{index}",
)
for index, center_y in enumerate(output_centers)
)
mode_monitors = (
bz.ModeMonitor(
center=(output_x, output_centers[1], 0.0),
size=(0.0, 3 * W1, domain_depth - 2 * pml_thickness),
freqs=freqs, mode_spec=mode_spec, name="mode_inner",
),
bz.ModeMonitor(
center=(output_x, output_centers[3], 0.0),
size=(0.0, 3 * W1, domain_depth - 2 * pml_thickness),
freqs=freqs, mode_spec=mode_spec, name="mode_outer",
),
)
field_monitor = bz.FieldMonitor(
center=(0.0, 0.0, 0.0), size=(domain_length, domain_width, 0.0),
freqs=[freq0], fields=("Ex", "Ey", "Ez", "Hx", "Hy", "Hz"), name="field",
)
monitors = (*flux_monitors, *mode_monitors, *((field_monitor,) if include_field else ()))
simulation = bz.Simulation(
domain=domain, design=design, sources=(source,), monitors=monitors,
boundaries=(bz.PML(thickness=pml_thickness, formulation="cpml"),),
grid_spec=grid_spec, run_time=(8 if test_mode else 5_000) / freq0,
setup_device="cpu" if test_mode else "default",
)
return simulation, {
"source": source, "output_x": output_x, "domain": domain,
"pml_thickness": pml_thickness, "mode_spec": mode_spec,
}3. Initial design and TM slot mode
The Tidy3D reference selects TM modes through mirror symmetry. Here the source and both mode monitors share an explicit polarization="tm" specification. Verify the fundamental mode's vertical Ez field, its peak in the GaN slot, and its even Ez parity about the horizontal and vertical mirrors before launching FDTD. The second returned mode is a higher TM-like mode; its index alone does not establish correspondence with the symmetry-filtered reference modes.
Plot only mode 0 for Ex, Ey, and Ez, with a separate magnitude scale for each component. BeamZ propagates along x, while the reference propagates along y: for magnitude comparisons, BeamZ Ex corresponds to reference Ey, BeamZ Ey corresponds to reference Ex, and Ez remains vertical in both. Compare field roles and component scales along with their shapes.
sim, setup = make_sim(L_MMI, W_MMI)
print(f"3D grid: {sim.grid.shape} | time steps: {sim.num_steps}")
print(f"Grid metrics: {sim.grid.metric_kind}; source frequency profiles: {len(setup['source'].profile_frequencies())}")
for axis in "xyz":
widths_nm = np.diff(sim.grid.axis_edges(axis)) / 1e-9
print(f"{axis} spacing: {widths_nm.min():.2f}–{widths_nm.max():.2f} nm")
fig, axes = sim.plot(
z=0.0, y=0.0, figsize=(12, 4.5), width_ratios=(1.2, 1.0),
source_markers=True, monitor_markers=True, show=False,
)
for axis in np.asarray(axes).flat:
axis.grid(False)
fig.suptitle("Initial 1×4 MMI: xy core plane and xz centerline", y=1.03)
plt.show()3D grid: (727, 239, 56) | time steps: 524905 Grid metrics: rectilinear; source frequency profiles: 5 x spacing: 27.78–27.79 nm y spacing: 25.00–66.21 nm z spacing: 25.00–65.94 nm
source_modes = setup["source"].solve_modes(sim, freqs=[freq0])
display(mode_data_to_dataframe(source_modes))
# Check the actual fundamental field, rather than trusting its mode index.
electric = np.asarray(source_modes.e_fields)[0, 0]
transverse_axes = source_modes.transverse_axes
edges = source_modes.grid_edges
centers = [
0.5 * (edge[:-1] + edge[1:]) - source_modes.center["xyz".index(axis)]
for axis, edge in zip(transverse_axes, edges, strict=True)
]
cell_area = np.multiply.outer(np.diff(edges[0]), np.diff(edges[1]))
ey_squared = np.abs(electric[1]) ** 2
ez_squared = np.abs(electric[2]) ** 2
tm_fraction = float(np.sum(ez_squared * cell_area) / np.sum((ey_squared + ez_squared) * cell_area))
peak_index = np.unravel_index(np.argmax(ez_squared), ez_squared.shape)
peak_z = float(centers[transverse_axes.index("z")][peak_index[transverse_axes.index("z")]])
ez = electric[2]
ez_parity = {}
for dimension, (axis, coordinates) in enumerate(zip(transverse_axes, centers, strict=True)):
reflected = np.apply_along_axis(
lambda row: np.interp(-coordinates, coordinates, row, left=0.0, right=0.0),
dimension, ez,
)
overlap = np.sum(np.conj(ez) * reflected * cell_area)
norm = np.sqrt(np.sum(np.abs(ez) ** 2 * cell_area) * np.sum(np.abs(reflected) ** 2 * cell_area))
ez_parity[axis] = float(np.real(overlap) / norm)
print(f"Fundamental TM check: vertical transverse-field fraction={tm_fraction:.4f}; peak z={peak_z / um:.4f} µm")
print("Even-Ez mirror overlaps:", {axis: round(value, 6) for axis, value in ez_parity.items()})
if not test_mode:
assert tm_fraction > 0.9, "The launched mode must be vertically polarized."
assert abs(peak_z) <= 0.5 * H_GAN, "The Ez peak must lie inside the GaN slot."
assert all(value > 0.95 for value in ez_parity.values()), "The fundamental Ez must have even mirror parity."
assert len(setup["source"].profile_frequencies()) == 5
fig, axes, neffs = plot_mode_field_components(
source_modes, field_names=("Ex", "Ey", "Ez"), mode_indices=(0,), val="abs",
f=freq0, figsize=(11.5, 3.8), cmap="magma", show=False,
)
for axis in np.asarray(axes).flat:
axis.grid(False)
axis.set_aspect("equal")
image = axis.collections[0]
image.set_clim(vmin=0.0)
fig.colorbar(image, ax=axis, label="Magnitude", fraction=0.046, pad=0.04)
fig.tight_layout()
plt.show()| f | mode_index | wavelength | n eff | k eff | loss (dB/cm) | mode area |
|---|---|---|---|---|---|---|
| 1.934145e+14 | 0 | 1.55 | 2.667952 | 0.0 | 0.0 | 0.196556 |
| 1 | 1.55 | 2.239197 | 0.0 | 0.0 | 0.323738 |
Fundamental TM check: vertical transverse-field fraction=0.9904; peak z=0.0375 µm
Even-Ez mirror overlaps: {'z': 0.999994, 'y': 1.0}
4. Run the initial design and inspect its response
The full run uses the direct CUDA backend. After the source has decayed, three consecutive checks must meet both the field-energy decay threshold (1e-5) and the relative DFT-monitor change threshold (1e-3). Checks cover all four output flux monitors and both mode monitors. The 1200-second budget includes compilation and chunked execution; an in-flight chunk and final result construction may exceed it. The reduced test uses a short fixed-duration JAX run.
from time import perf_counter
termination = None if test_mode else bz.AutoTermination(
max_wall_time=20 * 60,
field_decay=1e-5,
monitor_change=1e-3,
chunk_steps=2000,
consecutive_checks=3,
monitor_names=tuple(f"flux_{index}" for index in range(4)) + ("mode_inner", "mode_outer"),
)
wall_start = perf_counter()
initial_results = sim.run(
progress=not test_mode, backend=execution_backend, termination=termination,
)
print(f"Stopped: {initial_results.termination.reason if initial_results.termination else 'fixed duration'}")
if initial_results.termination is not None:
print(f"Measured {initial_results.termination.time / 1e-12:.3f} ps; converged: {initial_results.termination.converged}")
print(f"Wall time: {perf_counter() - wall_start:.1f} s (budget: 1200 s)")
spectrum_suffix = (
" — not yet converged"
if initial_results.termination is not None and not initial_results.termination.converged
else ""
)Simulation runtime: 1158.87 s GCUPS: 1.881 Stopped: wall_time_limit Measured 11.032 ps; converged: False Wall time: 1219.4 s (budget: 1200 s)
initial_results.plot_field("field", "Ez", frequency=freq0, val="real", cmap="RdBu", show_grid=False)
plt.title("Initial-design real Ez in the xy core plane")
plt.show()
def splitter_metrics(results):
launched = max(float(results.launched_power()), 1e-18)
output_power = np.stack([
np.maximum(np.asarray(results[f"flux_{index}"].flux, dtype=float), 0.0)
for index in range(4)
]) / launched
inner = 0.5 * (output_power[1] + output_power[2])
outer = 0.5 * (output_power[0] + output_power[3])
total = np.sum(output_power, axis=0)
excess_loss_db = -10 * np.log10(np.maximum(total, 1e-12))
imbalance = np.abs(inner - outer)
return output_power, inner, outer, excess_loss_db, imbalance
output_power, inner_power, outer_power, excess_loss_db, imbalance = splitter_metrics(initial_results)
center_index = int(np.argmin(np.abs(wavelengths - lambda0)))
fig, axes = plt.subplots(1, 2, figsize=(11, 4), constrained_layout=True)
for index, power in enumerate(output_power):
axes[0].plot(wavelengths / um, power, label=f"Output {index + 1}")
axes[0].set(xlabel="Wavelength (µm)", ylabel="Normalized output power", title="Transmission spectrum" + spectrum_suffix)
axes[0].legend(ncol=2, frameon=False)
axes[1].plot(wavelengths / um, excess_loss_db, color="#d81b60")
axes[1].axvline(lambda0 / um, color="0.35", ls="--", lw=1)
axes[1].set(xlabel="Wavelength (µm)", ylabel="Excess loss (dB)", title="Total excess loss" + spectrum_suffix)
for axis in axes:
axis.grid(False)
plt.show()
print(f"1550 nm: inner={inner_power[center_index]:.3f}, outer={outer_power[center_index]:.3f}, excess loss={excess_loss_db[center_index]:.3f} dB")
1550 nm: inner=0.237, outer=0.224, excess loss=0.354 dB
outer_mode_data = initial_results.mode("mode_outer")
forward_amps = np.asarray(outer_mode_data.amps.sel(direction="+"))
modal_fraction = np.abs(forward_amps) ** 2 / np.maximum(outer_power[:, None], 1e-12)
fig, ax = plt.subplots(figsize=(6, 3.8))
for mode_index in range(modal_fraction.shape[1]):
ax.plot(wavelengths / um, modal_fraction[:, mode_index], label=f"Mode {mode_index}")
ax.set(xlabel="Wavelength (µm)", ylabel="Outer-output modal power fraction", title="Mode composition" + spectrum_suffix)
ax.legend(frameon=False)
ax.grid(False)
plt.show()
# API opportunity: `results.mode(name).power_fraction(direction="+")`
# would eliminate this explicit amplitude-to-power normalization.
initial_results.plot_field("field", "E", frequency=freq0, val="abs^2", cmap="magma")
plt.show()