Skip to content
PhotonForge

MMI 3 dB Coupler with EME

Download notebook (.zip)

This example designs a 2x2 multimode interference (MMI) 3 dB coupler using the EMEModel, an eigenmode expansion (EME) solver.

MMI couplers are a textbook application of EME. Their operation relies on self-imaging: a field launched into a wide multimode region is periodically reproduced as single or multiple copies of the input as it propagates. EME is the natural method for this problem because it expands the field in the eigenmodes of each cross-section and propagates them analytically, which is far cheaper than time-stepping the whole volume with FDTD for a device that is long and dominated by propagation in a nearly uniform region.

The most valuable feature of EME for design is that the modes of a uniform section are solved only once, so the length of that section can be swept almost for free. We use exactly that capability to locate the 3 dB self-imaging length in a single simulation, and then visualize the resulting self-imaging field.

Technology

Section titled “Technology”

We use the built-in basic_technology, a silicon-on-insulator stack with a 220 nm device layer and a 500 nm wide single-mode strip waveguide (the "Strip" port specification).

import numpy as np
import photonforge as pf
import tidy3d as td
from matplotlib import pyplot as plt
from photonforge.live_viewer import LiveViewer


viewer = LiveViewer()


tech = pf.basic_technology()
pf.config.default_technology = tech


central_frequency = pf.C_0 / 1.55
LiveViewer started at http://127.0.0.1:46733. Call its 'stop()' method to free the port.

Parametric MMI component

Section titled “Parametric MMI component”

The component is a wide multimode region of length length and width mmi_width, fed by two input and two output access waveguides. Linear tapers widen the 500 nm access waveguides to taper_width at the junctions with the multimode region, which reduces excitation of unwanted modes and lowers the insertion loss. The geometry is generated with stencil.mmi and a cladding envelope is added so the Strip ports are detected correctly.

The access waveguides are placed at +/- port_separation / 2. Two points matter here:

  • The separation must be large enough that a port mode does not overlap the neighboring waveguide, otherwise the port mode is ill-defined. A center-to-center separation of 2 um (1 um edge-to-edge gap) is comfortable.
  • For a 2x2 coupler based on general interference, an input produces two self-images, one at the input transverse position and one mirrored about the center. Placing the outputs at the same +/- port_separation / 2 therefore captures both images.

We attach an EMEModel inside the function. The EME grid divides the device along the propagation direction into cells where modes are solved: a few cells resolve each taper, and a single cell covers the uniform multimode region, where one cell is exact because the propagation there is analytic. The accuracy and cost are set by num_modes and taper_cells; we justify the defaults of 20 and 8 in the convergence study near the end. An optional add_field_monitor flag attaches a small EMEFieldMonitor, restricted to the two low-order modes and the Ey component, so the propagated field can be recovered from a single-frequency run for visualization without producing a large dataset.

@pf.parametric_component
def mmi_coupler(
    *,
    length=56.5,
    mmi_width=4.0,
    port_length=5.0,
    taper_width=1.5,
    port_separation=2.0,
    num_modes=20,
    taper_cells=8,
    add_field_monitor=False,
):
    mmi = pf.Component("MMI_2x2")


    geometry = pf.stencil.mmi(
        length=length,
        width=mmi_width,
        num_ports=(2, 2),
        port_length=port_length,
        port_width=0.5,
        tapered_width=taper_width,
        port_separation=port_separation,
    )
    mmi.add("WG_CORE", *geometry)
    mmi.add("WG_CLAD", pf.envelope(mmi, offset=1.0, trim_x_min=True, trim_x_max=True))
    mmi.add_port(mmi.detect_ports([tech.ports["Strip"]]))


    # EME cells: taper_cells per taper, a single cell for the uniform region.
    x_min = mmi.bounds()[0][0]
    mode_spec = td.EMEModeSpec(num_modes=num_modes)
    eme_grid = td.EMECompositeGrid(
        subgrids=[
            td.EMEUniformGrid(num_cells=taper_cells, mode_spec=mode_spec),
            td.EMEUniformGrid(num_cells=1, mode_spec=mode_spec),
            td.EMEUniformGrid(num_cells=taper_cells, mode_spec=mode_spec),
        ],
        subgrid_boundaries=[x_min + port_length, x_min + port_length + length],
    )


    monitors = []
    if add_field_monitor:
        # Keep the field data small: only the two low-order (even/odd) modes and
        # the Ey component are needed to reconstruct a single-port TE field.
        monitors = [
            td.EMEFieldMonitor(
                center=(0, 0, 0.11),
                size=(td.inf, td.inf, 0),
                name="field",
                num_modes=2,
                fields=["Ey"],
            )
        ]
    mmi.add_model(pf.EMEModel(eme_grid_spec=eme_grid, monitors=monitors), "EME")


    return mmi


mmi = mmi_coupler()
viewer(mmi)

The four ports are detected as P0, P1 (inputs, left) and P2, P3 (outputs, right). With a single input at P0, P2 is the bar output (same side as the input) and P3 is the cross output.

Finding the 3 dB length with an EME length sweep

Section titled “Finding the 3 dB length with an EME length sweep”

The length of the uniform multimode region sets the splitting ratio. A first estimate comes from the beat length of the two lowest modes of the multimode region, L_pi = pi / (beta_0 - beta_1). For this 4 um wide region L_pi is about 39 um, and the two-fold (3 dB) self-image of general interference forms near 3 * L_pi / 2, about 59 um. This is only a starting point, so we sweep the length to locate the true 3 dB point.

Instead of rebuilding and re-simulating the device at every length, we use the native EME length sweep through the simulation_updates argument of the model. A tidy3d.EMELengthSweep scales the length of selected EME cells without re-solving the modes: only the analytic propagation phase changes. We scale only the single central cell, so every length in the sweep reuses the same modes and the whole sweep runs as one simulation.

The scale_factors array has one row per swept length and one column per EME cell. All columns are 1 except the central-cell column, which carries length / base_length.

base_length = 56.0
swept_lengths = np.linspace(20.0, 75.0, 56)


mmi_sweep = mmi_coupler(length=base_length)
eme_grid = mmi_sweep.models["EME"].eme_grid_spec


subgrids = eme_grid.subgrids
num_cells = sum(sub.num_cells for sub in subgrids)
center_cell = subgrids[0].num_cells  # index of the uniform-region cell


# Scale only the central (uniform) cell, leaving the taper cells fixed.
scale_factors = np.ones((swept_lengths.size, num_cells))
scale_factors[:, center_cell] = swept_lengths / base_length


# A dedicated model carries the length sweep so the component's own model,
# reused later for characterization, is left untouched.
sweep_model = pf.EMEModel(
    eme_grid_spec=eme_grid,
    simulation_updates={"sweep_spec": td.EMELengthSweep(scale_factors=scale_factors.tolist())},
)
eme_data = sweep_model.simulation_data(mmi_sweep, [central_frequency], show_progress=False)
Uploading task 'EME…'
Uploading task 'Mode-Stripwaveguide…'
Uploading task 'Mode-Stripwaveguide…'
Uploading task 'Mode-Stripwaveguide…'
Uploading task 'Mode-Stripwaveguide…'
Starting task 'EME': https://tidy3d.simulation.cloud/workbench?taskId=eme-541056cf-bf95-4023-b1b4-fa851be39b64
Starting task 'Mode-Stripwaveguide': https://tidy3d.simulation.cloud/workbench?taskId=mo-5ce23979-98e0-4b4b-806c-06633eef4b62
Starting task 'Mode-Stripwaveguide': https://tidy3d.simulation.cloud/workbench?taskId=mo-fb66be0b-f371-42d8-b728-7b43c1526fc8
Starting task 'Mode-Stripwaveguide': https://tidy3d.simulation.cloud/workbench?taskId=mo-f8bc6aa2-a6ba-4f7d-83fb-eb0cc7612b48
Downloading data from 'EME'…
Starting task 'Mode-Stripwaveguide': https://tidy3d.simulation.cloud/workbench?taskId=mo-140c6f9f-7823-48de-adb2-2cb7d6b94e5a
Downloading data from 'Mode-Stripwaveguide'…
Downloading data from 'Mode-Stripwaveguide'…
Downloading data from 'Mode-Stripwaveguide'…
Downloading data from 'Mode-Stripwaveguide'…

The swept scattering matrix eme_data.smatrix is expressed in the basis of the EME port modes. Because the two access waveguides are well separated, the port modes of each facet are the symmetric (even) and antisymmetric (odd) supermodes of the waveguide pair, returned as mode indices 0 and 1. A single-waveguide excitation is a balanced combination of these supermodes,

P0=12 (even+odd),P2=12 (even+odd),P3=12 (even−odd),P_0 = \tfrac{1}{\sqrt 2}\,(\text{even} + \text{odd}), \qquad P_2 = \tfrac{1}{\sqrt 2}\,(\text{even} + \text{odd}), \qquad P_3 = \tfrac{1}{\sqrt 2}\,(\text{even} - \text{odd}),

so the bar and cross transmissions for a P0 input follow from the four supermode-to-supermode terms of S21. This even/odd decomposition is the standard way to read physical-port results from a 2x2 MMI EME simulation (see the Tidy3D community example Cascaded 2x2 MMI array); it reproduces what the high-level s_matrix call computes internally. We wrap it in a small helper so the same step serves both the length sweep here and the mode-convergence sweep later.

def bar_cross(s21):
    """Physical bar and cross power for a P0 input, from the EME supermode
    S21 block. The last two axes of `s21` are (mode_out, mode_in)."""
    ee, eo = s21[..., 0, 0], s21[..., 0, 1]
    oe, oo = s21[..., 1, 0], s21[..., 1, 1]
    bar = np.abs(0.5 * (ee + eo + oe + oo)) ** 2
    cross = np.abs(0.5 * (ee + eo - oe - oo)) ** 2
    return bar, cross


t_bar, t_cross = bar_cross(eme_data.smatrix.S21.isel(f=0).values)


# 3 dB length: the most balanced split among low-loss, long-device solutions.
candidate = (swept_lengths > 45.0) & (t_bar + t_cross > 0.85)
best = np.argmin(np.where(candidate, np.abs(t_bar - t_cross), np.inf))
optimal_length = swept_lengths[best]
print(f"3 dB length: {optimal_length:.1f} um")
print(f"  bar = {t_bar[best]:.3f}, cross = {t_cross[best]:.3f}, "
      f"total = {t_bar[best] + t_cross[best]:.3f}")
3 dB length: 57.0 um
  bar = 0.486, cross = 0.487, total = 0.972
fig, ax = plt.subplots(figsize=(7, 4), tight_layout=True)
ax.plot(swept_lengths, t_bar, label="bar (P0 to P2)")
ax.plot(swept_lengths, t_cross, label="cross (P0 to P3)")
ax.axvline(optimal_length, color="k", ls="--", lw=1)
ax.set_xlabel("Multimode region length (um)")
ax.set_ylabel("Transmitted power")
ax.set_title("EME length sweep (single simulation)")
ax.legend()
ax.grid(True, alpha=0.3)
plt.show()

The bar and cross curves cross near the predicted length, where the input power splits equally between the two outputs. Around 28 um the device acts instead as a cross-over (the single mirrored image), and the 3 dB point sits at roughly twice that length, as expected for two-fold imaging.

Convergence

Section titled “Convergence”

The cost and accuracy of the model are controlled by the number of modes kept in each cell and the number of cells resolving the tapers. We verify both at the design length. The mode count is swept in a single simulation with tidy3d.EMEModeSweep, which reuses one mode solve and repeats only the propagation with fewer modes. The taper-cell count changes the mode-solving grid, so it needs a separate build per value.

# Mode-count convergence in a single simulation with EMEModeSweep.
mode_counts = [4, 8, 12, 16, 20, 24]
conv = mmi_coupler(length=optimal_length, num_modes=max(mode_counts))
mode_sweep_model = pf.EMEModel(
    eme_grid_spec=conv.models["EME"].eme_grid_spec,
    simulation_updates={"sweep_spec": td.EMEModeSweep(num_modes=mode_counts)},
)
mode_s21 = mode_sweep_model.simulation_data(
    conv, [central_frequency], show_progress=False
).smatrix.S21.isel(f=0).values
mode_bar, mode_cross = bar_cross(mode_s21)


# Taper-cell convergence: one build per value (the mode grid changes).
cell_counts = [2, 4, 8, 12]
cell_total = []
for n in cell_counts:
    s = mmi_coupler(length=optimal_length, taper_cells=n).s_matrix(
        [central_frequency], model_kwargs={"inputs": ["P0"]}
    )
    cell_total.append(np.abs(s["P0@0", "P2@0"][0]) ** 2 + np.abs(s["P0@0", "P3@0"][0]) ** 2)


fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 3.2), tight_layout=True)
ax1.plot(mode_counts, mode_bar + mode_cross, "o-")
ax1.set_xlabel("Number of modes")
ax1.set_ylabel("Total transmission")
ax1.set_title("Mode convergence (EMEModeSweep)")
ax1.grid(True, alpha=0.3)
ax2.plot(cell_counts, cell_total, "s-")
ax2.set_xlabel("Cells per taper")
ax2.set_ylabel("Total transmission")
ax2.set_title("Taper-cell convergence")
ax2.grid(True, alpha=0.3)
plt.show()
Uploading task 'EME…'
Uploading task 'Mode-Stripwaveguide…'
Uploading task 'Mode-Stripwaveguide…'
Uploading task 'Mode-Stripwaveguide…'
Uploading task 'Mode-Stripwaveguide…'
Starting task 'Mode-Stripwaveguide': https://tidy3d.simulation.cloud/workbench?taskId=mo-53ff32f9-d405-461e-80b7-a553381247c2
Starting task 'Mode-Stripwaveguide': https://tidy3d.simulation.cloud/workbench?taskId=mo-007119c8-3371-4c4e-aa00-3af7eaf3c25a
Starting task 'Mode-Stripwaveguide': https://tidy3d.simulation.cloud/workbench?taskId=mo-e5d97fc8-8998-47ef-9b7c-7a5fce27d54d
Starting task 'EME': https://tidy3d.simulation.cloud/workbench?taskId=eme-64ada359-9cd2-4582-9cbb-32ea078562f8
Starting task 'Mode-Stripwaveguide': https://tidy3d.simulation.cloud/workbench?taskId=mo-f9d29103-4177-485d-8f44-34d00a0d5e7c
Downloading data from 'Mode-Stripwaveguide'…
Downloading data from 'Mode-Stripwaveguide'…
Downloading data from 'Mode-Stripwaveguide'…
Downloading data from 'EME'…
Downloading data from 'Mode-Stripwaveguide'…
Uploading task 'EME…'
Starting task 'EME': https://tidy3d.simulation.cloud/workbench?taskId=eme-c620b3dc-87e9-4704-bc67-a07d9c5291e8
Downloading data from 'EME'…
Uploading task 'EME…'
Starting task 'EME': https://tidy3d.simulation.cloud/workbench?taskId=eme-cd73516f-c94a-4471-a685-20d8ef4ce7ab
Downloading data from 'EME'…
Uploading task 'EME…'
Starting task 'EME': https://tidy3d.simulation.cloud/workbench?taskId=eme-57f1e892-5a41-48da-aadb-a674a25383d1
Downloading data from 'EME'…
Uploading task 'EME…'
Starting task 'EME': https://tidy3d.simulation.cloud/workbench?taskId=eme-cf588e12-f6e0-4be1-bf26-7f124e200d8b
Downloading data from 'EME'…

Both curves flatten by num_modes = 20 and taper_cells = 8, confirming the defaults are converged while keeping the simulation small.

Visualizing the self-imaging

Section titled “Visualizing the self-imaging”

We recover the field from a small dedicated run at the center wavelength: the coupler is rebuilt with add_field_monitor=True, and its monitor stores only the two low-order modes and the Ey component, so a single frequency keeps the dataset small. The field for a single P0 input follows from the same even/odd combination used for the scattering matrix.

field_coupler = mmi_coupler(length=optimal_length, add_field_monitor=True)
field_data = field_coupler.models["EME"].simulation_data(
    field_coupler, [central_frequency], show_progress=False
)["field"]


ey = field_data.Ey.isel(eme_port_index=0, sweep_index=0, z=0, f=0)
p0_ey = (ey.isel(mode_index=0).values + ey.isel(mode_index=1).values) / np.sqrt(2)
intensity = np.nan_to_num(np.abs(p0_ey) ** 2)


x = field_data.Ey.x.values
y = field_data.Ey.y.values


fig, ax = plt.subplots()


ax.pcolormesh(x, y, intensity.T / intensity.max(), cmap="magma_r", shading="auto")


# Scale the y-direction by a factor of 8
ax.set_aspect(5)


ax.set_ylim(-2.5, 2.5)
ax.set_xlabel("x (um)")
ax.set_ylabel("y (um)")
ax.set_title("Self-imaging field intensity for a P0 input")
plt.show()
Uploading task 'EME…'
Starting task 'EME': https://tidy3d.simulation.cloud/workbench?taskId=eme-2d66ac98-6373-43aa-8141-cc9af3266724
Downloading data from 'EME'…

Light launched into the lower input refocuses into two equal images at the output facet, which is the self-imaging that produces the balanced 3 dB split.

At the 3 dB length the input power splits almost equally between the two outputs with low excess loss, as the length sweep reports and the self-imaging field confirms. As a 2x2 MMI, the coupler also imparts the characteristic 90 degree phase between its bar and cross outputs, the signature of a 3 dB hybrid. The design was cross-checked against a full 3D FDTD simulation (Tidy3DModel), which reproduces the same balanced split and quadrature phase.

Reference

  1. Soldano, Lucas B., and Erik C. M. Pennings. “Optical multi-mode interference devices based on self-imaging: principles and applications.” Journal of Lightwave Technology 1995 13 (4), 615-627, doi: 10.1109/50.372474.
  2. Tidy3D community example, “Cascaded 2x2 MMI array: EME simulation and S-parameter extraction,” flexcompute.com/tidy3d/community/notebooks/MMIEME.