Cascaded MZI Filter

This example shows the design of a cascaded Mach-Zehnder interferometer (MZI) filter based on [1]. The 4-stage filter is constructed from directional couplers with specific ratios of power coupling and preset phase delays between the arms for each stage.

We start by designing the 5 directional couplers using a parametric component function and finding the required coupling lengths to reach each power coupling ratio. Then we design the 4 arm sections from the effective and group indices fro the waveguide profile we’re using. Finally, we test the complete filter and investigate how changes in the fabrication might impact the filter response.

References

1. 19. Dwivedi, P. De Heyn, P. Absil, J. Van Campenhout and W. Bogaerts, “Coarse wavelength division multiplexer on silicon-on-insulator for 100 GbE,” 2015 IEEE 12th International Conference on Group IV Photonics (GFP), Vancouver, BC, Canada, 2015, pp. 9-10, doi: 10.1109/Group4.2015.7305928.

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

We use the SiEPIC OpenEBL technology to design this filter. Let’s take a look at the layers and port specifications provided by this technology.

import siepic_forge as siepic

tech = siepic.ebeam()

print("Ports:")
for k, v in sorted(tech.ports.items()):
    print(f"- {k}: {v.width} μm, {v.num_modes} mode(s)")

tech.layers

07:24:54 -03 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.                                       

Ports:
- MM_TE_1550_2000: 6.0 μm, 10 mode(s)
- MM_TE_1550_3000: 6.0 μm, 15 mode(s)
- Rib_TE_1310_350: 2.35 μm, 1 mode(s)
- Rib_TE_1550_500: 2.5 μm, 1 mode(s)
- Slot_TE_1550_500: 3.0 μm, 1 mode(s)
- TE-TM_1550_450: 2.2 μm, 2 mode(s)
- TE_1310_350: 1.5 μm, 1 mode(s)
- TE_1310_410: 1.5 μm, 1 mode(s)
- TE_1550_500: 1.5 μm, 1 mode(s)
- TM_1310_350: 1.5 μm, 1 mode(s)
- TM_1550_500: 1.5 μm, 1 mode(s)
- eskid_TE_1550: 2.0 μm, 1 mode(s)

Name	Layer	Description	Color	Pattern
Si	(1, 0)	SiEPIC - Waveguide	#ff80a818	\\
PinRec	(1, 10)	SiEPIC	#ff80a818	xx
PinRecM	(1, 11)	SiEPIC	#80000018	+
Si Slab	(2, 0)	Dedicated Run Layers - Device… … Layer Partial Etch	#c080ff18	/
Direct Metal	(5, 0)	Dedicated Run Layers	#80a8ff18	||
Oxide open to BOX	(6, 0)	Dedicated Run Layers	#ff000018	-
Text	(10, 0)	Text-Not Fabricated	#00000018	hollow
M1_heater	(11, 0)	TiW Heater	#0000ff18	\\
M2_router	(12, 0)	TiW/Au Routing Bilayer	#ffbf0018	//
M_Open	(13, 0)	Bond Pad Open	#80005718	\\
Si n	(20, 0)	Dedicated Run Layers	#afff8018	-
Si p	(21, 0)	Dedicated Run Layers	#ffd9df18	=
Si n+	(22, 0)	Dedicated Run Layers	#ff800018	x
Si p+	(23, 0)	Dedicated Run Layers	#ddff0018	xx
Si n++	(24, 0)	Dedicated Run Layers	#00ffff18	+
Si p++	(25, 0)	Dedicated Run Layers	#00800018	++
ANT Reserved	(31, 0)	SiEPIC/ANT Reserved	#9580ff18	/
ANT Reserved 1	(33, 0)	ANT Reserved	#9580ff18	/
Via to silicon	(40, 0)	Dedicated Run Layers	#0000ff18	.
DevRec	(68, 0)	SiEPIC	#00800018	.
FbrTgt	(81, 0)	SiEPIC/Dedicated Run Layers	#80808018	++
ANT Reserved 2	(102, 0)	ANT Reserved	#9580ff18	/
ANT Reserved 3	(110, 0)	ANT Reserved	#9580ff18	/
Custom Dicing	(189, 0)		#00000018	hollow
SEM Imaging	(200, 0)		#ff000018	x
Deep Trench	(201, 0)		#00ff0018	.
Deep Trench Handling Exclusion	(202, 0)		#00760018	:
Thermal Isolation Trenches	(203, 0)		#00800018	\
Laser Integration Shelf	(205, 0)	Dedicated Run Layers	#69ff0518	xx
Floor Plan-Not Fabricated	(290, 0)		#c080ff18	hollow
Error: device layer width is… … less than design rule	(301, 0)	DRC Errors	#80005718	=
Error: device layer spacing is… … less than design rule	(301, 1)	DRC Errors	#80005718	-
Warning: polygons/paths on… … PinRec layer (1/10) will NOT be fabricated	(301, 2)	DRC Errors	#80005718	||
Error: direct metal width is… … less than 5 microns	(305, 0)	DRC Errors	#80808018	++
Error: direct metal spacing is… … less than 10 microns	(305, 1)	DRC Errors	#80808018	+
Error: TiW width is less than 3… … microns	(311, 0)	DRC Errors	#ffa08018	//
Error: TiW spacing is less than… … 3 microns	(311, 1)	DRC Errors	#ffa08018	/
Error: Al width is less than… … design rule	(312, 0)	DRC Errors	#00ffff18	|
Error: Al spacing is less than… … design rule	(312, 1)	DRC Errors	#00ffff18	//
Error: Spacing between TiW and… … Al is less than 5 microns	(312, 3)	DRC Errors	#00ffff18	\\
Error: Oxide window width is… … less than 10 microns	(313, 0)	DRC Errors	#01ff6b18	||
Error: Oxide window spacing is… … less than 10 microns	(313, 1)	DRC Errors	#01ff6b18	|
Error: Oxide window is not… … placed over Al	(313, 2)	DRC Errors	#01ff6b18	//
Standard Design Area	(350, 0)	DRC Errors	#ddff0018	\
Error: Features outside design… … area. Verify design size and centering.	(350, 1)	DRC Errors	#ddff0018	:
Error: Dicing lane width is… … less than 100 microns	(389, 0)	DRC Errors	#ff00ff18	++
Error: Spacing between dicing… … lane and devices is less than 50 microns	(389, 1)	DRC Errors	#ff00ff18	+
Error: SEM width is less than… … 500 nm	(400, 0)	DRC Errors	#ff9d9d18	x
Deep Trench Design Area	(401, 0)	DRC Errors	#80a8ff18	xx
Error: Metal, SEM, or handling… … region overlap with deep trenches. Verify design centering	(401, 1)	DRC Errors	#80a8ff18	x
Warning: Silicon features… … outside deep trench design area. Verify accuracy before submission	(401, 2)	DRC Errors	#80a8ff18	=
Error: Spacing between metal… … and deep trench is less than 30 microns	(401, 3)	DRC Errors	#80a8ff18	-
Error: Deep trench width is… … less than 260 microns	(401, 4)	DRC Errors	#80a8ff18	||
Error: Deep trench handling… … area missing. Please add handling area of size shown by polygons	(402, 0)	DRC Errors	#ff000018	+
Error: Features inside deep… … trench handling area	(402, 1)	DRC Errors	#ff000018	xx
Error: Thermal isolation width… … is less than design rule	(403, 0)	DRC Errors	#50008018	++
Error: Thermal isolation… … spacing is less than design rule	(403, 1)	DRC Errors	#50008018	+
Error: Spacing between thermal… … isolation and metal is less than design rule	(403, 2)	DRC Errors	#50008018	xx
Error: Thermal isolation and… … device layer overlap, or spacing is less than design rule	(403, 3)	DRC Errors	#50008018	x
Dream Photonics Black Box-Not… … Fabricated	(998, 0)		#00000018	hollow
Errors	(999, 0)	SiEPIC	#0000ff18	||

We define our wavelength range of interest and select a port specification appropriate for a single-mode TE device:

pf.config.default_technology = tech
wavelengths = np.linspace(1.53, 1.63, 11)

port_spec = "TE_1550_500"

# Get the dimensions of our main waveguide geometry from the port specification
wg_width, _ = tech.ports[port_spec].path_profile_for("Si")
clad_width = (tech.ports[port_spec].width - wg_width) / 2

The filter specifications are chosen based on the main reference [1].

# Free spectral range
fsr = 0.01

# Defined in the reference work
power_ratios = [0.5, 0.13, 0.12, 0.5, 0.25]

We also define default arguments for parametric components:

pf.config.default_kwargs = {
    "port_spec": port_spec,
    "radius": 10,
    "euler_fraction": 0.5,
    "tidy3d_model_kwargs": {"verbose": False},
    "waveguide_model_kwargs": {"verbose": False},
}

Coupler Design

We use the parametric ring_coupler for the basic coupler design. The only parameter of note here is the coupling distance, which we define as the waveguide core width plus the coupling gap.

offset = clad_width / 2
coupling_gap = 0.12

coupler = pf.parametric.s_bend_coupler(
    coupling_distance=wg_width + coupling_gap,
    coupling_length=1,
    s_bend_offset=offset,
    model=pf.Tidy3DModel(
        port_symmetries=[
            ("P1", "P0", "P3", "P2"),  # x-axis reflection
            ("P2", "P3", "P0", "P1"),  # y-axis reflection
            ("P3", "P2", "P1", "P0"),  # inversion symmetry
        ]
    ),
)

coupler

# Test the defined port symmetries (assert if test fails)
assert coupler.active_model.test_port_symmetries(
    coupler, frequencies=[pf.C_0 / 1.55], grid_spec=8
)

Uploading task 'P0@0…'
Uploading task 'P1@0…'
Uploading task 'P2@0…'
Uploading task 'P3@0…'
Starting task 'P3@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-e1da14b0-ab4a-466f-b8a0-02567d6e7972
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-5e3231b3-60d9-4b38-9784-9428f8ce0c5d
Downloading data from 'P3@0'…
Downloading data from 'P0@0'…
Starting task 'P1@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-b921d522-4c16-49e0-9824-a39a47e3a146
Starting task 'P2@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-9e830808-4602-441d-a8a6-f383df882756
Downloading data from 'P1@0'…
Downloading data from 'P2@0'…
Progress: 100%
Progress: 100%

We want to test our coupler to make sure the losses and reflections are low. Because we have included a model in our parametric component, we can look at its scattering parameters in the wavelength range of interest directly with the plot_s_matrix utility.

_ = pf.plot_s_matrix(coupler.s_matrix(td.C_0 / wavelengths), input_ports=["P0"])

Uploading task 'P0@0…'
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-96676827-2667-41f0-ba4f-89863ca43ad5
Downloading data from 'P0@0'…
Progress: 100%

Losses are hard to gauge directly in those plots. Since we are also interested in the power coupling ratio, we can define a helper function to calculate those figures directly for our couplers:

def get_coupler_metrics(coupling_length):
    coupler.update(coupling_length=coupling_length)
    s_matrix = coupler.s_matrix(td.C_0 / wavelengths, model_kwargs={"inputs": ["P0@0"]})
    thru = np.abs(s_matrix[("P0@0", "P2@0")]) ** 2
    drop = np.abs(s_matrix[("P0@0", "P3@0")]) ** 2
    trans = drop + thru
    return drop / trans, 1 - trans


_, ax = plt.subplots(1, 1, figsize=(5, 3.5), tight_layout=True)
ratio, loss = get_coupler_metrics(1)
ax.plot(wavelengths, ratio, label="Power ratio")
ax.plot(wavelengths, loss, label="Loss")
ax.set(ylim=(0, None), xlabel="λ (μm)")
_ = ax.legend()

Progress: 100%

To design each coupler based on a desired power ratio, we will use a simple sweep of coupling lengths and linear interpolation to find the coupling length required for the each ratio. We will sweep a few values of coupling length until the ratio starts to decrease or we have surpassed the maximal required power ratio.

Note that, depending on the waveguide materials, coupling gap, and S bend shape, we might not be able to reach all required power ratios. In that case, those other parameters have to be changed, or the coupler redesigned.

sweep_lengths = [0.0]
sweep_ratios = []
sweep_losses = []
w = wavelengths.size // 2
while True:
    print(f"Testing coupling length = {sweep_lengths[-1]}…")
    ratio, loss = get_coupler_metrics(coupling_length=sweep_lengths[-1])
    sweep_ratios.append(ratio)
    sweep_losses.append(loss)
    if (
        len(sweep_ratios) > 1
        and (sweep_ratios[-1][w] < sweep_ratios[-2][w])
        or sweep_ratios[-1][w] > max(*power_ratios)
    ):
        break
    sweep_lengths.append(sweep_lengths[-1] + 1)

sweep_lengths = np.array(sweep_lengths)
sweep_losses = np.array(sweep_losses)
sweep_ratios = np.array(sweep_ratios)

_, ax = plt.subplots(1, 2, figsize=(10, 3.5), tight_layout=True)

for i in [round(x) for x in np.linspace(0, wavelengths.size - 1, 5)]:
    ax[0].plot(sweep_lengths, sweep_ratios[:, i], label=f"{wavelengths[i]:.2f} µm")

for length, loss in zip(sweep_lengths, sweep_losses):
    ax[1].plot(wavelengths, loss, label=f"{length:.1f} µm")

ax[0].set(ylim=(0, None), xlabel="Coupler length (μm)", ylabel="Power ratio")
ax[0].legend()
ax[1].set(ylim=(0, None), xlabel="Wavelength (μm)", ylabel="Loss")
_ = ax[1].legend()

Testing coupling length = 0.0…
Uploading task 'P0@0…'Progress: 0% \
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-85713b9e-d418-4b77-8d96-3dc1377a8428
Downloading data from 'P0@0'…
Progress: 100%
Testing coupling length = 1.0…
Starting…

07:25:49 -03 Loading simulation from local cache. View cached task using web UI
             at
             'https://tidy3d.simulation.cloud/workbench?taskId=fdve-96676827-266
             7-41f0-ba4f-89863ca43ad5'.

Progress: 100%
Testing coupling length = 2.0…
Uploading task 'P0@0…'
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-e9c348b6-7d94-40d2-88bd-4dd68a05411a
Downloading data from 'P0@0'…
Testing coupling length = 3.0…
Progress: 100%
Uploading task 'P0@0…'
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-c5c69418-64a2-4653-afc6-a4000f13e6c2
Downloading data from 'P0@0'…
Progress: 100%
Testing coupling length = 4.0…
Uploading task 'P0@0…'Progress: 0% \
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-cf5329ec-353a-451e-9072-46f48de4bc25
Downloading data from 'P0@0'…
Progress: 100%
Testing coupling length = 5.0…
Uploading task 'P0@0…'
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-22daed02-4233-4687-9efe-6a1bf4c2ba9e
Downloading data from 'P0@0'…
Progress: 100%
Testing coupling length = 6.0…
Uploading task 'P0@0…'Progress: 0% \
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-dae1a837-eab2-428a-b7ad-b6e9a29a1002
Downloading data from 'P0@0'…
Progress: 100%

We use the gathered ratios and coupling lengths to back-calculate the length from the ratio using linear interpolation. Note that the result is rounded to 10 nm to be easier to inspect (and because any rounding error will not have significant impact in the final device).

def get_coupler_length(desired_ratio):
    pr = sweep_ratios[:, wavelengths.size // 2]
    if any(np.diff(pr) < 0):
        m = np.where(np.diff(pr) < 0)[0][0] + 1
    else:
        m = len(sweep_lengths)
    return np.round(np.interp(desired_ratio, pr[:m], sweep_lengths[:m]), decimals=2)


get_coupler_length(power_ratios)

array([5.47, 1.28, 1.13, 5.47, 2.84])

At this point we are ready to create the couplers corresponding to each power ratio required by the filter. We compute and plot their responses to make sure they are within specification.

dcs = {
    pr: coupler.copy(True).update(coupling_length=get_coupler_length(pr))
    for pr in set(power_ratios)
}
dcs[power_ratios[0]]

_, ax = plt.subplots(1, 2, figsize=(10, 3.5), tight_layout=True)

for pr in sorted(set(power_ratios)):
    ratio, loss = get_coupler_metrics(get_coupler_length(pr))
    ax[0].plot(wavelengths, ratio, label=str(pr))
    ax[1].plot(wavelengths, loss)
    ax[0].axhline(pr, c="tab:gray", lw=0.5)

ax[0].axvline(wavelengths.mean(), c="tab:gray", lw=0.5)
ax[0].set(ylim=(0, None), xlabel="λ (μm)", ylabel="Power ratio")
ax[1].set(ylim=(0, None), xlabel="λ (μm)", ylabel="Loss")
_ = ax[0].legend()

Uploading task 'P0@0…'Progress: 0% \
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-019ebb84-594e-4f15-9430-ee194279a1e0
Downloading data from 'P0@0'…
Progress: 100%
Uploading task 'P0@0…'
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-77c4f947-a813-46f5-8e89-d78eac2592e8
Downloading data from 'P0@0'…
Progress: 100%
Uploading task 'P0@0…'
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-2cb503b4-b9ed-4229-8f32-a467c7f66156
Downloading data from 'P0@0'…
Progress: 100%
Uploading task 'P0@0…'
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-3dfd0a5d-81f9-4ef9-bd9a-a59cd422dcd6
Downloading data from 'P0@0'…
Progress: 100%

Phase Delays

The relation between the arm lengths, the phase delays, and free spectral range (FSR) of the filter require the knowledge of the effective and group indices of the waveguide mode we are using. That information can be quickly obtained through the port_modes function.

mode_solver = pf.port_modes(
    tech.ports[port_spec],
    [td.C_0 / wavelengths[wavelengths.size // 2]],
    mesh_refinement=40,
    group_index=True,
)
mode_solver.data.to_dataframe()

Uploading task 'Mode-ModeSolver…'
Starting task 'Mode-ModeSolver': https://tidy3d.simulation.cloud/workbench?taskId=mo-9b6fbb2a-2920-4268-b92d-5d5da5436ca4
Downloading data from 'Mode-ModeSolver'…
Progress: 100%

		wavelength	n eff	k eff	loss (dB/cm)	TE (Ey) fraction	wg TE fraction	wg TM fraction	mode area	group index	dispersion (ps/(nm km))
f	mode_index
1.897421e+14	0	1.58	2.40887	0.0	0.0	0.981923	0.75657	0.81593	0.200152	4.190646	452.25128

Following the main reference [1], we compute the lengths required for a 180° phase shift and for the FSR. With those, the path differences for each stage are calculated, again, rounding the results to reasonable values.

dl = wavelengths[wavelengths.size // 2] ** 2 / (
    fsr * float(mode_solver.data.n_group.isel(f=0, mode_index=0, drop=True))
)
l_pi = wavelengths[wavelengths.size // 2] / (
    2 * float(mode_solver.data.n_eff.isel(f=0, mode_index=0, drop=True))
)

path_differences = np.round([dl, 2 * dl, l_pi - 2 * dl, -2 * dl], decimals=2).tolist()
path_differences

[59.57, 119.14, -118.81, -119.14]

We create both arms of an MZI stage as a single component. This option makes it easier to implement both arms bending to the same side with some separation between the external and internal arms.

Similar to what we did for the coupler, we create a 90° bend that will be used for all arms, as well as a section of a straight waveguide used as a separator.

Because the arms are, ideally, independent of each other, we can use a circuit model to compute its S parameters based on the models for the bend and straight sections. The straight and bend sections can be modeled with a semi-analytical waveguide model, which is accurate for isolated waveguides and efficient to compute independent of the waveguide length. Note that it only works well for the bend because it has a large enough radius to make radiation losses negligible.

separation = (
    dcs[power_ratios[0]].ports["P1"].center[1]
    - dcs[power_ratios[0]].ports["P0"].center[1]
)

separator = pf.parametric.straight(length=separation, name="SEPARATOR")

bend = pf.parametric.bend(name="BEND")


@pf.parametric_component
def mzi_arms(*, path_difference):
    delta = pf.parametric.straight(
        length=path_difference / 2, name=f"DELTA_{path_difference}"
    )

    c = pf.Component(f"ARMS_{path_difference}")

    arm = c.add_reference(separator)
    c.add_port(arm["P0"], "P0")
    arm = c.add_reference(bend).connect("P0", arm["P1"])
    arm = c.add_reference(bend).connect("P1", arm["P1"])
    arm = c.add_reference(bend).connect("P1", arm["P0"])
    arm = c.add_reference(bend).connect("P0", arm["P0"])
    arm = c.add_reference(separator).connect("P0", arm["P1"])
    c.add_port(arm["P1"], "P2")

    arm = c.add_reference(bend).translate((0, separation))
    c.add_port(arm["P0"], "P1")
    arm = c.add_reference(delta).connect("P0", arm["P1"])
    arm = c.add_reference(bend).connect("P1", arm["P1"])
    arm = c.add_reference(separator).connect("P0", arm["P0"])
    arm = c.add_reference(separator).connect("P0", arm["P1"])
    arm = c.add_reference(bend).connect("P1", arm["P1"])
    arm = c.add_reference(delta).connect("P0", arm["P0"])
    arm = c.add_reference(bend).connect("P0", arm["P1"])
    c.add_port(arm["P1"], "P3")

    # Set verbose=True to see the upload and download progress for each simulation
    c.add_model(pf.CircuitModel(verbose=False), "Circuit")
    return c


mzi_arms(path_difference=10)

Our parametric arm component only works for positive path differences, so we generate all required arms based on the absolute value of the differences. Note that 2 of the sections differ only in sign (second and fourth stages), so the same sub-component can be used for both, but mirrored on the stage with negative sign.

arms = {dl: mzi_arms(path_difference=dl) for dl in {abs(pd) for pd in path_differences}}
arms[abs(path_differences[0])]

Complete Filter

Once we have all sub-components that compose the filter, it is easy to assemble them together. We only have to be careful with required rotation of phase delay arms with negative signs.

The complete filter is quite large (we can check its bounds to verify it is larger than 200 × 160 μm²). Although possible, simulating the whole geometry in Tidy3D requires a very refined mesh to guarantee the phase delays are accurately accumulated in each stage. Another option is to use the circuit model to combine the independent responses from all sub-components. We use this option to take advantage of the semi-analytical waveguide models and smaller simulation sizes for the bend and couplers.

cascaded_mzi = pf.Component("MZI_FILTER")

coupler_ref = cascaded_mzi.add_reference(dcs[power_ratios[0]])
cascaded_mzi.add_port(coupler_ref["P0"])
cascaded_mzi.add_port(coupler_ref["P1"])

for pr, dl in zip(power_ratios[1:], path_differences):
    if dl > 0:
        p0, p2 = "P0", "P2"
    else:
        p0, p2 = "P3", "P1"
    arm_ref = cascaded_mzi.add_reference(arms[abs(dl)]).connect(p0, coupler_ref["P2"])
    coupler_ref = cascaded_mzi.add_reference(dcs[pr]).connect("P0", arm_ref[p2])

cascaded_mzi.add_port(coupler_ref["P2"])
cascaded_mzi.add_port(coupler_ref["P3"])

cascaded_mzi.add_model(pf.CircuitModel(verbose=False), "Circuit")
cascaded_mzi.write_gds()

We compute the final scattering parameters using a finer wavelength spacing than the one used for coupler design because we want to capture the finer features in the filter response.

dense_wavelengths = np.linspace(1.53, 1.63, 401)

s_filter = cascaded_mzi.s_matrix(td.C_0 / dense_wavelengths)

Uploading task 'P0@0…'
Uploading task 'Mode-StripTE1550nmw500nm…'
Uploading task 'Mode-StripTE1550nmw500nm…'
Uploading task 'P0@0…'
Uploading task 'P0@0…'
Uploading task 'P0@0…'
Starting task 'Mode-StripTE1550nmw500nm': https://tidy3d.simulation.cloud/workbench?taskId=mo-fe03172a-980e-40c1-bfab-8ead23b0fe3c
Downloading data from 'Mode-StripTE1550nmw500nm'…
Starting task 'Mode-StripTE1550nmw500nm': https://tidy3d.simulation.cloud/workbench?taskId=mo-a82f8731-de9b-42ff-8b6a-836823629f0d
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-6bc58ea4-e827-485c-9f4b-bb56281efe86
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-0631172c-47c9-4914-807d-50feeb665545
Downloading data from 'Mode-StripTE1550nmw500nm'…
Downloading data from 'P0@0'…
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-710a04c2-9345-403d-aa83-90873c904ce7
Downloading data from 'P0@0'…
Downloading data from 'P0@0'…
Starting task 'P0@0': https://tidy3d.simulation.cloud/workbench?taskId=fdve-0b9b0785-3cd8-4ab8-b559-43527772463d
Downloading data from 'P0@0'…
Progress: 100%

fig, ax = plt.subplots(1, 1, figsize=(12, 4))

ax.plot(
    dense_wavelengths, 20 * np.log10(np.abs(s_filter[("P0@0", "P2@0")])), label="Thru"
)
ax.plot(
    dense_wavelengths, 20 * np.log10(np.abs(s_filter[("P0@0", "P3@0")])), label="Cross"
)

ax.set(ylim=(-40, None), xlabel="λ (µm)", ylabel="Transmission (dB)")
_ = ax.legend()

Monte Carlo Analysis

We can perform a simple Monte Carlo analysis by introducing a length error in the delay lines. This analysis is very fast because the lengths affect only semi-analytical waveguide models, so we will use it as an example (it is also somewhat related to geometrical imperfections leading to phase errors, so it is a valid analysis). In practice, it might be more interesting to explore the effects of variations in the technology, such as a material layer thickness or some mask dilation, as long as the technology exposes them as parameters.

Technologies in PhotonForge can provide preset random variables for their parameters based on measured data, which can greatly accelerate the setup of the Monte Carlo analysis. The technology we’re using includes 1 such variable: it describes the thickness of the device silicon layer as a normal distribution with mean of 220 nm and standard deviation of approximately 4 nm.

tech.random_variables

[RandomVariable('si_thickness', **{'value': 0.22, 'stdev': 0.0037166666666666667})]

For this example, we have to create appropriate a list of RandomVariable for the waveguide models. The original length parameter will be used as the mean value for the normal distribution, with a standard deviation of 0.05%. We can get the original length value directly from the parametric keyword arguments defined in the component.

random_variables = []
for component in cascaded_mzi.dependencies():
    if component.name.startswith("ARMS"):
        length = component.parametric_kwargs["path_difference"]
        random_variables.append(
            pf.monte_carlo.RandomVariable(
                "path_difference", component, value=length, stdev=5e-4 * length
            )
        )

random_variables

[RandomVariable('path_difference', **{'value': 118.81, 'stdev': 0.059405}),
 RandomVariable('path_difference', **{'value': 119.14, 'stdev': 0.059570000000000005}),
 RandomVariable('path_difference', **{'value': 59.57, 'stdev': 0.029785000000000002})]

The s_matrix function within the monte_carlo sub-module can be used to automate the analysis and return the S parameters for a number of random samples and corner cases once we define the random variables to be used. We collect 12 random cases in this example:

_, results = pf.monte_carlo.s_matrix(
    cascaded_mzi,
    td.C_0 / dense_wavelengths,
    *random_variables,
    random_samples=12,
    random_seed=0,
)

Starting sample 1 of 12…
Starting sample 2 of 12…
Starting sample 3 of 12…
Starting sample 4 of 12…
Starting sample 5 of 12…
Starting sample 6 of 12…
Starting sample 7 of 12…
Starting sample 8 of 12…
Starting sample 9 of 12…
Starting sample 10 of 12…
Starting sample 11 of 12…
Starting sample 12 of 12…
Sample 1 done.
Sample 2 done.
Sample 3 done.
Sample 4 done.
Sample 5 done.
Sample 6 done.
Sample 7 done.
Sample 8 done.
Sample 9 done.
Sample 10 done.
Sample 11 done.
Sample 12 done.
All samples done!

The parameter values for each sample and the corresponding S matrices are returned in a list:

print("Samples:", len(results))
print("Variable values for the first sample:", results[0][:-1])

Samples: 12
Variable values for the first sample: (59.57446634723455, 119.21801089942575, 118.81526395129751)

fig, ax = plt.subplots(1, 1, figsize=(12, 4))

for *_, s_matrix in results:
    ax.plot(
        dense_wavelengths,
        20 * np.log10(np.abs(s_matrix[("P0@0", "P2@0")])),
        alpha=0.4,
        color="tab:blue",
    )

_ = ax.set(ylim=(-40, None), xlabel="λ (µm)", ylabel="Transmission (dB)")

Variations by Reference

One issue with the previous analysis is apparent from the difference in the numbers of interferometric arms and random variables. That’s because the random variables are applied to component parameters, and 2 of the arms are references to the same component, as can be seen in the reference list:

cascaded_mzi.references[3], cascaded_mzi.references[7]

(Reference(component="ARMS_119.14", origin=(61.602, 0), rotation=0, scaling=1, x_reflection=False, repetition=Repetition(columns=1, rows=1, spacing=(0, 0))),
 Reference(component="ARMS_119.14", origin=(207.534, 1.12), rotation=-180, scaling=1, x_reflection=False, repetition=Repetition(columns=1, rows=1, spacing=(0, 0))))

cascaded_mzi.references[3].component is cascaded_mzi.references[7].component

True

That means that the variations applied to the second and fourth arms are identical, because they come from a single variable and are applied to the component itself. In physical terms, it is as if the variations in those arms are perfectly correlated, which is not the analysis that we intended.

Instead, we want to modify the component independently in each reference to it. We can achieve that when using a Circuit model through by-reference updates via the updates argument in its start function. This argument is used by the Monte Carlo runner when we define random variables with reference specifications.

Each definition is a tuple with the reference specification and a dictionary with random variables that should be used to update the component of the specified reference, its technology or its model. The specification is another tuple that describes the “path” to the reference, and is detailed in the Circuit model documentation for the start function.

Let’s create 4 random variables, one for each arm using their references. Our specification will be the component name (we know they are direct references of the main component, so there’s no need to dive into the dependency tree) and the index for that component, identifying the number of occurrences for each component.

reference_variables = []
for pd in path_differences:
    dl = abs(pd)
    ref_spec = (f"ARMS_{dl}", 0)

    # make sure we increment the index for the second reference to the component used twice
    if any(ref_spec == r for r, _ in reference_variables):
        ref_spec = (f"ARMS_{dl}", 1)

    updates = {
        "component_updates": [
            pf.monte_carlo.RandomVariable("path_difference", value=dl, stdev=5e-4 * dl)
        ]
    }
    reference_variables.append((ref_spec, updates))

reference_variables

[(('ARMS_59.57', 0),
  {'component_updates': [RandomVariable('path_difference', **{'value': 59.57, 'stdev': 0.029785000000000002})]}),
 (('ARMS_119.14', 0),
  {'component_updates': [RandomVariable('path_difference', **{'value': 119.14, 'stdev': 0.059570000000000005})]}),
 (('ARMS_118.81', 0),
  {'component_updates': [RandomVariable('path_difference', **{'value': 118.81, 'stdev': 0.059405})]}),
 (('ARMS_119.14', 1),
  {'component_updates': [RandomVariable('path_difference', **{'value': 119.14, 'stdev': 0.059570000000000005})]})]

Note that we could also apply the updates to the waveguide models of each “DELAY” reference, if, for example, the arms were not created as parametric components. For that, the following variables could be used:

reference_variables_2 = []
for pd in path_differences:
    dl = abs(pd)
    ref_spec = (f"ARMS_{dl}", 0, "DELTA_.*")

    # Make sure we increment the index for the second reference to the component used twice
    if any(ref_spec == r for r, _ in reference_variables_2):
        ref_spec = (f"ARMS_{dl}", 1, "DELTA_.*")

    # Note that each delay section is half the length of the full path difference
    updates = {
        "model_updates": [
            pf.monte_carlo.RandomVariable("length", value=dl / 2, stdev=5e-4 * dl / 2)
        ]
    }
    reference_variables_2.append((ref_spec, updates))

reference_variables_2

[(('ARMS_59.57', 0, 'DELTA_.*'),
  {'model_updates': [RandomVariable('length', **{'value': 29.785, 'stdev': 0.014892500000000001})]}),
 (('ARMS_119.14', 0, 'DELTA_.*'),
  {'model_updates': [RandomVariable('length', **{'value': 59.57, 'stdev': 0.029785000000000002})]}),
 (('ARMS_118.81', 0, 'DELTA_.*'),
  {'model_updates': [RandomVariable('length', **{'value': 59.405, 'stdev': 0.0297025})]}),
 (('ARMS_119.14', 1, 'DELTA_.*'),
  {'model_updates': [RandomVariable('length', **{'value': 59.57, 'stdev': 0.029785000000000002})]})]

We will use the component updates in this case, but feel free to test with the model ones. Because of the predefined random seed, the results will be identical!

_, results = pf.monte_carlo.s_matrix(
    cascaded_mzi,
    td.C_0 / dense_wavelengths,
    *reference_variables,
    random_samples=12,
    random_seed=0,
)

Starting sample 1 of 12…
Starting sample 2 of 12…
Starting sample 3 of 12…
Starting sample 4 of 12…
Starting sample 5 of 12…
Starting sample 6 of 12…
Starting sample 7 of 12…
Starting sample 8 of 12…
Starting sample 9 of 12…
Starting sample 10 of 12…
Starting sample 11 of 12…
Starting sample 12 of 12…
Sample 1 done.
Sample 2 done.
Sample 3 done.
Sample 4 done.
Sample 5 done.
Sample 6 done.
Sample 7 done.
Sample 8 done.
Sample 9 done.
Sample 10 done.
Sample 11 done.
Sample 12 done.
All samples done!

fig, ax = plt.subplots(1, 1, figsize=(12, 4))

for *_, s_matrix in results:
    ax.plot(
        dense_wavelengths,
        20 * np.log10(np.abs(s_matrix[("P0@0", "P2@0")])),
        alpha=0.4,
        color="tab:blue",
    )

_ = ax.set(ylim=(-40, None), xlabel="λ (µm)", ylabel="Transmission (dB)")