Quantum Chip

This notebook is an example for full-chip layout and routing using PhotonForge. The layout is inspired by the two-qubit processor presented by Qiang et al. in [1]. The exact design is not available in the reference, but we can adapt it based on the figures in the paper.

This re-implementation is based on the open source LNOI400 PDK from Luxtelligence.

References

1. Qiang, X., et al. “Large-scale silicon quantum photonics implementing arbitrary two-qubit processing,” Nature Photon 12, 534–539 (2018), doi: 10.1038/s41566-018-0236-y.

import luxtelligence_lnoi400_forge as lxt
import numpy as np
import photonforge as pf
import tidy3d as td
from photonforge import live_viewer

viewer = live_viewer.LiveViewer()

Starting live viewer at http://localhost:5001

We start by defining a default argument for parametric components, including the default technology, that will be used throughout the design.

tech = lxt.lnoi400(include_substrate=False, include_top_opening=False)
pf.config.default_technology = tech

port_spec = tech.ports["RWG1000"]

pf.config.default_kwargs = {
    "port_spec": port_spec,
    "radius": 75,
    "euler_fraction": 0.5,
}

Individual components

A few small components will be used repeatedly in the layout, so we can start by designing each of those. Each component can be simulated and optimized as demonstrated in other examples. To keep this notebook more focused on layout, we will assume that optimization has been done separately and no simulations will be run here.

The waveguide crossing and spiral can be derived directly from the internal parametric library:

cross = pf.parametric.crossing(arm_length=40, added_width=2)
viewer(cross)

spiral = pf.parametric.rectangular_spiral(full_length=12_000)

print(f"Spiral length: {pf.route_length(spiral) * 1e-3:.6f} mm")
viewer(spiral)

Spiral length: 11.999995 mm

The 2×2 MMI comes from the LNOI400 PDK:

mmi = lxt.component.mmi2x2()
viewer(mmi)

We will also need a termination element to avoid extraneous reflections in unused ports. We design the termination as a slightly bend waveguide that is tapered down to radiate light efficiently with minimal reflections. Again, we assume the proper optimizations have been already made externally.

The implementation is in the form of a custom parametric component, to facilitate optimization and Monte Carlo studies.

@pf.parametric_component
def termination(*, length, offset_factor=0.3, sign=1):
    c = pf.Component()

    for layer, path in port_spec.get_paths((0, 0)):
        w = None
        l = length

        if layer == tech.layers["LN_RIDGE"].layer:
            w = pf.config.tolerance
            l -= 0.5 * port_spec.width

        c.add(
            layer, path.bezier([(0.5 * l, 0), (l, sign * offset_factor * l)], width=w)
        )

    c.add_port(pf.Port((0, 0), 0, port_spec))

    # We use an analitical termination model to remove all reflections, but optimization
    # of the design could be done with a Tidy3D model.
    c.add_model(pf.TerminationModel())
    return c


terminations = [
    termination(length=2 * pf.config.default_kwargs["radius"], sign=-1),
    termination(length=2 * pf.config.default_kwargs["radius"], sign=1),
]

viewer(terminations[0])

Phase-controlled MZI

Mach-Zehnder interferometers (MZI) with phase shift control are key building blocks in the quantum processor. Before diving into the MZI itself, we can create the thermal-controlled phase-shifter as a separate component, including its own optical ports, electrical terminals, and model.

@pf.parametric_component
def phase_shifter(
    *,
    length=100,
    heater_width=3,
    heater_pad_size=25,
    heater_pad_distance=30,
):
    c = pf.Component()
    ref = c.add_reference(pf.parametric.straight(length=length))
    c.add_port([ref["P0"], ref["P1"]])

    heater = pf.Rectangle(center=(0.5 * length, 0), size=(length, heater_width))
    c.add("HT", heater)

    pad = pf.Rectangle(size=(heater_pad_size, heater_pad_size))
    pad0 = pf.Terminal("HT", pad.copy().translate((-heater_pad_distance, 0)))
    pad1 = pf.Terminal("HT", pad.copy().translate((length + heater_pad_distance, 0)))
    c.add_terminal([pad0, pad1])

    connection = pf.Rectangle(size=(heater_width, heater_width))
    conn0 = pf.Terminal("HT", connection)
    conn1 = pf.Terminal("HT", connection.copy().translate((length, 0)))
    c.add(
        pf.parametric.route_taper(terminal1=pad0, terminal2=conn0),
        pf.parametric.route_taper(terminal1=pad1, terminal2=conn1),
    )

    # The waveguide model can be used to emulate heating by artificially changing
    # its length parameter later.
    c.add_model(pf.WaveguideModel())

    return c


ps = phase_shifter()
viewer(ps)

We will export a GDSII layout for the phase shifter because it will be used in the Tunable MZI example. Note that a better option for re-use would be to store either a phf file or storing the parametric function as a custom PDK that can be loaded later.

ps.write_gds()

The MZI parametric component is designed to offer enough flexibility to cover all use cases that appear in the quantum processor, including unbalanced arms, phase-shift control, chained stages, etc. The choice of making a single function to cover all of those cases allows us to avoid code duplication and maintain consistency throughout the chip. However, the resulting function becomes more complex with every new argument, so a balance must be reached in a case-by-case basis.

@pf.parametric_component
def mzi(
    *,
    added_length=0,
    stages=1,
    upper_arm="wg",
    lower_arm="ps",
    output_mmi=True,
):
    """
    MZI component.

    Args:
        added_length: Waveguide length added to the upper arm (negative values
          add to the lower arm).
        stages: Number of MZI stages.
        upper_arm: Upper arm specification can be ``'wg'``, ``'ps'``, or
          ``None``, for a bare waveguide, a phase-controlled arm, or no
          arm, respectively.
        lower_arm: Lower arm specification, analogous to ``uper_arm``.
        output_mmi: If ``False``, the output MMI in the last stage is
          suppressed.
    """
    c = pf.Component()

    # Phase-shifter arm
    ps = phase_shifter()
    ps_len = ps.parametric_kwargs["length"]

    # Bare waveguide arm
    wg = pf.parametric.straight(length=ps_len)

    # Minimal arm offsets and distance w.r.t. the MMIs
    offset = max(2 * ps.parametric_kwargs["heater_pad_size"], port_spec.width)

    if added_length == 0:
        # Balanced case: use S bends for routing
        route = pf.parametric.route_s_bend
        distance = 6 * offset
        # Offsets for lower and upper arms
        offsets = [-offset, offset]
    else:
        # Increase offset to guarantee that no S bends are used in the route,
        # which, in turn, guarantees that the added length is exact.
        distance = 2 * pf.config.default_kwargs["radius"]
        offset = max(offset, distance + pf.config.tolerance)
        route = pf.parametric.route
        # Offsets for lower and upper arms
        offsets = [min(0, added_length / 2) - offset, max(0, added_length / 2) + offset]

    # Calculate position to align MMIs at the x axis
    x0 = mmi.ports["P0"].center[0]
    y0 = 0.5 * (mmi.ports["P1"].center[1] + mmi.ports["P0"].center[1])

    # For multi-stage, calculate spacing
    mmi_len = mmi.ports["P2"].center[0] - mmi.ports["P0"].center[0]
    spacing = (mmi_len + ps_len + 2 * distance, 0)
    mmi_cols = (stages + 1) if output_mmi else stages

    # Add all MMIs
    mmi_ref = pf.Reference(mmi, origin=(-x0, -y0), columns=mmi_cols, spacing=spacing)
    c.add(mmi_ref)

    # Add ports. If no output MMI, add output ports from arms later.
    c.add_port([mmi_ref.get_ports("P0")[0], mmi_ref.get_ports("P1")[0]])
    if output_mmi:
        c.add_port([mmi_ref.get_ports("P2")[-1], mmi_ref.get_ports("P3")[-1]])

    # Arm start position
    x_ps = mmi_ref.get_ports("P2")[0].center[0] + distance

    # Add each arm to all stages
    for i, arm in enumerate([lower_arm, upper_arm]):
        arm = {"wg": wg, "ps": ps}.get(arm)
        if arm is None:
            # Use terminations on both MMIs
            mmi_port = mmi_ref.get_ports(f"P{2 + i}")[0]
            term_ref = pf.Reference(terminations[i], columns=stages, spacing=spacing)
            term_ref.connect("P0", mmi_port)
            c.add(term_ref)

            # There's only 1 MMI for a single stage without output MMI
            if mmi_cols > 1:
                mmi_port = mmi_ref.get_ports(f"P{i}")[1]
                term_ref = pf.Reference(
                    terminations[1 - i], columns=mmi_cols - 1, spacing=spacing
                )
                term_ref.connect("P0", mmi_port)
                c.add(term_ref)
        else:
            arm = pf.Reference(arm, (x_ps, offsets[i]), columns=stages, spacing=spacing)
            rt = pf.Reference(
                route(port1=(mmi_ref, f"P{2 + i}", 0), port2=(arm, "P0", 0)),
                columns=stages,
                spacing=spacing,
            )
            c.add(arm, rt)

            # Route to next stage only when there is a second MMI
            if mmi_cols > 1:
                rt = pf.Reference(
                    route(port1=(mmi_ref, f"P{i}", 1), port2=(arm, "P1", 0)),
                    columns=mmi_cols - 1,
                    spacing=spacing,
                )
                c.add(rt)

            # Add output port if no output MMI is used
            if not output_mmi:
                c.add_port(arm.get_ports("P1")[-1])

            # Add terminals, if using a phase-shifter
            if len(arm.component.terminals) > 0:
                arm_terminals = arm.get_terminals()
                c.add_terminal(
                    [
                        arm_terminals[name][index]
                        for index in range(stages)
                        for name in ("T0", "T1")
                    ]
                )

    c.add_model(pf.CircuitModel())

    return c


viewer(mzi())

viewer(mzi(stages=3))

viewer(mzi(added_length=100))

viewer(mzi(lower_arm=None))

viewer(mzi(output_mmi=False))

One structure that appears repeatedly in the design is a two-stage, phase-controlled MZI without output MMI. We create it to be used later whenever needed:

coupling_mzi = mzi(stages=2, output_mmi=False)
viewer(coupling_mzi)

Chip layout

Being a complex design, we split the full device into sub-components that can be re-used and combined into the full layout later. Each is implemented as a parametric component to enable advanced analyses later, if necessary.

@pf.parametric_component
def input_section(*, mzi_unbalance=-200):
    c = pf.Component()

    mzi_in = c.add_reference(coupling_mzi)
    c.add_terminal([mzi_in["T0"], mzi_in["T1"], mzi_in["T2"], mzi_in["T3"]])

    s_ref0 = c.add_reference(spiral).mirror().connect("P0", mzi_in["P2"])
    s_ref1 = c.add_reference(spiral).connect("P0", mzi_in["P3"])

    unbalanced = mzi(
        added_length=mzi_unbalance, output_mmi=False, upper_arm="ps", lower_arm="wg"
    )

    # Connect references to align ports
    mzi_u0 = c.add_reference(unbalanced).connect("P0", s_ref0["P1"])
    mzi_u1 = c.add_reference(unbalanced).connect("P0", s_ref1["P1"])

    # We need room for terminations and to avoid any overlaps between the MZIs
    gap = terminations[0].parametric_kwargs["length"] + port_spec.width
    if mzi_u0.y_max > 0 or mzi_u1.y_min < 0:
        gap += 2 * pf.config.default_radius - port_spec.width
        dy = max(mzi_u0.y_max, -mzi_u1.y_min)
        mzi_u0.y_mid -= dy
        mzi_u1.y_mid += dy
    mzi_u0.x_min += gap
    mzi_u1.x_min += gap

    # Add routes and terminations
    c.add(
        pf.parametric.route(port1=(s_ref0, "P1"), port2=(mzi_u0, "P0")),
        pf.parametric.route(port1=(s_ref1, "P1"), port2=(mzi_u1, "P0")),
        pf.Reference(terminations[0]).connect("P0", mzi_u0["P1"]),
        pf.Reference(terminations[0]).connect("P0", mzi_u1["P1"]),
    )

    c.add_terminal([mzi_u0["T0"], mzi_u0["T1"], mzi_u1["T0"], mzi_u1["T1"]])

    single = mzi(upper_arm=None, lower_arm="wg")
    distance = 2 * pf.config.default_kwargs["radius"]
    mzi_out0 = c.add_reference(single).translate(
        (mzi_u0.x_max + distance, mzi_u0.origin[1])
    )
    mzi_out1 = c.add_reference(single).translate(
        (mzi_u1.x_max + distance, mzi_u1.origin[1])
    )

    c.add(pf.parametric.route(port1=(mzi_u0, "P2"), port2=(mzi_out0, "P0")))
    c.add(pf.parametric.route(port1=(mzi_u0, "P3"), port2=(mzi_out0, "P1")))
    c.add(pf.parametric.route(port1=(mzi_u1, "P2"), port2=(mzi_out1, "P0")))
    c.add(pf.parametric.route(port1=(mzi_u1, "P3"), port2=(mzi_out1, "P1")))

    c.add_port(
        [
            mzi_in["P0"],
            mzi_in["P1"],
            mzi_out0["P2"],
            mzi_out0["P3"],
            mzi_out1["P2"],
            mzi_out1["P3"],
        ]
    )
    c.add_model(pf.CircuitModel())

    return c


in_section = input_section()
viewer(in_section)

@pf.parametric_component
def dual_mzi_array(*, stages=6):
    c = pf.Component()

    # The output stages are different for each row
    mzi_out = [
        mzi(lower_arm="wg", upper_arm="ps", output_mmi=False),
        mzi(lower_arm="ps", upper_arm="wg", output_mmi=False),
    ]

    mzi_offset = mzi_out[0].ports["P3"].center[1] - mzi_out[0].ports["P2"].center[1]
    cross_size = cross.ports["P3"].center[1] - cross.ports["P1"].center[1]
    offset = 2 * pf.config.default_kwargs["radius"] + cross_size + 2 * mzi_offset

    mzi_row = pf.Reference(
        mzi(stages=stages - 1, output_mmi=False), rows=2, spacing=(0, offset)
    )
    c.add(mzi_row)
    ports = mzi_row.get_ports()
    c.add_port([ports[name][index] for index in range(2) for name in ["P0", "P1"]])
    terminals = mzi_row.get_terminals()

    x_out = (
        mzi_row.get_ports("P0")[0].center[0]
        + (stages - 1) * mzi_row.component.references[0].spacing[0]
    )

    for i, component in enumerate(mzi_out):
        ref = pf.Reference(component, (x_out, offset * i))
        c.add(
            ref,
            pf.parametric.route_s_bend(port1=(mzi_row, "P2", i), port2=(ref, "P0")),
            pf.parametric.route_s_bend(port1=(mzi_row, "P3", i), port2=(ref, "P1")),
        )
        c.add_port([ref["P2"], ref["P3"]])
        c.add_terminal(
            [terminals[f"T{j * 2 + k}"][i] for j in range(stages - 1) for k in range(2)]
        )
        c.add_terminal([ref["T0"], ref["T1"]])

    c.add_model(pf.CircuitModel())

    return c


mzi_array = dual_mzi_array()
viewer(mzi_array)

@pf.parametric_component
def output_section():
    c = pf.Component()

    distance = 2 * pf.config.default_kwargs["radius"]
    cross_size = cross.ports["P3"].center[1] - cross.ports["P1"].center[1]
    offset = mzi_array.ports["P2"].center[1] - mzi_array.ports["P0"].center[1]

    array = pf.Reference(mzi_array, rows=2, spacing=(0, 2 * offset))
    array_ports = array.get_ports()
    array_terminals = array.get_terminals()
    c.add(array)
    c.add_port(
        [
            array_ports[name][index]
            for index in range(2)
            for name in ("P0", "P1", "P2", "P3")
        ]
    )
    c.add_terminal(
        [
            array_terminals[f"T{j}"][index]
            for index in range(2)
            for j in range(len(array_terminals))
        ]
    )

    # Add space for crossings between arrays and output MMIs
    x_mmi = (
        array_ports["P4"][0].center[0]
        + 4 * distance
        + 3 * cross_size
        - mmi.ports["P0"].center[0]
    )
    y_mmi = array_ports["P0"][0].center[1] + offset / 2 - mmi.ports["P0"].center[1]
    mmi4 = pf.Reference(mmi, (x_mmi, y_mmi), rows=4, spacing=(0, offset))
    mmi4_ports = mmi4.get_ports()
    c.add(mmi4)

    c.add(pf.parametric.route(port1=(mmi4, "P0", 0), port2=(array, "P4", 0)))
    c.add(pf.parametric.route(port2=(mmi4, "P1", 3), port1=(array, "P7", 1)))

    x_cross = array_ports["P4"][0].center[0] + distance + cross_size / 2
    y_cross = (array_ports["P5"][0].center[1] + array_ports["P6"][0].center[1]) / 2
    cross0 = pf.Reference(cross, (x_cross, y_cross), rows=3, spacing=(0, offset))
    c.add(cross0)

    c.add(
        pf.parametric.route(port1=(cross0, "P1", i * 2), port2=(array, "P5", i))
        for i in range(2)
    )
    c.add(
        pf.parametric.route(port1=(cross0, "P0", i * 2), port2=(array, "P6", i))
        for i in range(2)
    )
    c.add(
        pf.parametric.route(port1=(cross0, "P1", 1), port2=(array, "P7", 0)),
        pf.parametric.route(port1=(cross0, "P0", 1), port2=(array, "P4", 1)),
        pf.parametric.route(port2=(mmi4, "P1", 0), port1=(cross0, "P2", 0)),
        pf.parametric.route(
            port1=(mmi4, "P0", 3),
            port2=(cross0, "P3", 2),
            waypoints=[mmi4_ports["P0"][3].center - pf.config.default_kwargs["radius"]],
        ),
    )

    x_cross += distance + cross_size
    y_cross += offset
    cross1 = pf.Reference(cross, (x_cross, y_cross), rows=2, spacing=(0, offset))
    c.add(cross1)

    c.add(
        pf.parametric.route_s_bend(port1=(cross0, "P3", i), port2=(cross1, "P1", i))
        for i in range(2)
    )
    c.add(
        pf.parametric.route(port1=(cross0, "P2", i + 1), port2=(cross1, "P0", i))
        for i in range(2)
    )
    c.add(
        pf.parametric.route(port1=(cross1, "P2", 0), port2=(mmi4, "P0", 1)),
        pf.parametric.route(port1=(mmi4, "P1", 2), port2=(cross1, "P3", 1)),
    )

    x_cross += distance + cross_size
    y_cross += offset / 2
    cross2 = pf.Reference(cross, (x_cross, y_cross))
    c.add(cross2)

    c.add(
        pf.parametric.route(port1=(cross2, "P0"), port2=(cross1, "P3", 0)),
        pf.parametric.route(port1=(cross2, "P3"), port2=(cross1, "P2", 1)),
        pf.parametric.route(port1=(cross2, "P2"), port2=(mmi4, "P0", 2)),
        pf.parametric.route(port1=(cross2, "P1"), port2=(mmi4, "P1", 1)),
    )

    # Second stage of MMIs
    x_mmi = mmi4_ports["P2"][0].center[0] + distance - mmi.ports["P0"].center[0]
    y_mmi = mmi4_ports["P0"][0].center[1] + offset / 2 - mmi.ports["P0"].center[1]
    mmi2 = pf.Reference(mmi, (x_mmi, y_mmi), rows=2, spacing=(0, 2 * offset))
    mmi2_ports = mmi2.get_ports()
    c.add(mmi2)

    c.add(
        pf.parametric.route(port1=(mmi4, "P2", i * 2 + 1), port2=(mmi2, "P1", i))
        for i in range(2)
    )
    c.add(
        pf.parametric.route(port1=(mmi4, "P3", i * 2), port2=(mmi2, "P0", i))
        for i in range(2)
    )

    for i in range(2):
        c.add_reference(terminations[0]).connect("P0", mmi4_ports["P2"][i * 2])
        c.add_reference(terminations[1]).connect("P0", mmi4_ports["P3"][i * 2 + 1])

    # Output MZI
    x_mzi_out = mmi2_ports["P2"][0].center[0] + distance + coupling_mzi["P2"].center[0]
    y_mzi_out = mmi2_ports["P0"][0].center[1] + offset
    mzi_out = pf.Reference(coupling_mzi, (x_mzi_out, y_mzi_out), 180, x_reflection=True)
    c.add(mzi_out)
    c.add_port([mzi_out["P0"], mzi_out["P1"]])
    c.add_terminal([mzi_out["T3"], mzi_out["T2"], mzi_out["T1"], mzi_out["T0"]])

    c.add(
        pf.parametric.route(port1=(mmi2, "P2", 1), port2=(mzi_out, "P3")),
        pf.parametric.route(port1=(mmi2, "P3", 0), port2=(mzi_out, "P2")),
    )

    c.add_reference(terminations[0]).connect("P0", mmi2_ports["P2"][0])
    c.add_reference(terminations[1]).connect("P0", mmi2_ports["P3"][1])

    c.add_model(pf.CircuitModel())

    return c


out_section = output_section()
viewer(out_section)

Final layout

We combine each sub-component into the final layout and route all terminals to side pads.

@pf.parametric_component
def quantum_processor(*, pad_size=100):
    c = pf.Component()

    # Input coupler
    mzi_in = c.add_reference(coupling_mzi)
    c.add_port([mzi_in["P0"], mzi_in["P1"]])

    # Use the terminal size for spacing metal wires
    (x_min, _), (x_max, _) = coupling_mzi.terminals["T0"].bounds()
    w = 2 * (x_max - x_min)

    (_, y_min), (_, y_max) = out_section.bounds()
    output_size = y_max - y_min

    cross_size = cross.ports["P3"].center[1] - cross.ports["P1"].center[1]

    distance = 2 * pf.config.default_kwargs["radius"]
    offset = (
        output_size
        + coupling_mzi.ports["P3"].center[1]
        - coupling_mzi.ports["P2"].center[1]
    )

    # Two input sections routed to the input coupler
    ref_in = pf.Reference(
        in_section,
        (mzi_in["P2"].center[0] + distance, -offset / 2),
        rows=2,
        spacing=(0, offset),
    )
    c.add(ref_in)

    c.add(
        pf.parametric.route(port1=(mzi_in, "P2"), port2=(ref_in, "P1", 0)),
        pf.parametric.route(port1=(mzi_in, "P3"), port2=(ref_in, "P0", 1)),
    )

    c.add_reference(terminations[1]).connect("P0", ref_in["P0"][0])
    c.add_reference(terminations[0]).connect("P0", ref_in["P1"][1])

    # Two output sections
    x_out = ref_in["P2"][0].center[0] + 4 * distance + 3 * cross_size
    y_out = -(out_section["P7"].center[1] + out_section["P0"].center[1] + offset) / 2
    ref_out = pf.Reference(out_section, (x_out, y_out), rows=2, spacing=(0, offset))
    c.add(ref_out)
    c.add_port([ref_out["P8"][0], ref_out["P9"][0], ref_out["P8"][1], ref_out["P9"][1]])

    c.add(
        pf.parametric.route(port1=(ref_in, "P2", 0), port2=(ref_out, "P1", 0)),
        pf.parametric.route(port1=(ref_in, "P5", 1), port2=(ref_out, "P7", 1)),
    )

    for i in range(8):
        c.add_reference(terminations[1]).connect("P0", ref_out[f"P{i // 2 * 2}"][i % 2])

    # Routing between input and output sections using 3 crossing stages

    x_cross = ref_in["P2"][0].center[0] + distance + cross_size / 2
    y_cross = -offset / 2
    cross0 = pf.Reference(cross, (x_cross, y_cross), rows=3, spacing=(0, offset / 2))
    c.add(cross0)

    c.add(
        pf.parametric.route(port1=(cross0, "P1", 0), port2=(ref_in, "P3", 0)),
        pf.parametric.route(port1=(cross0, "P1", 1), port2=(ref_in, "P4", 0)),
        pf.parametric.route(port1=(cross0, "P1", 2), port2=(ref_in, "P3", 1)),
        pf.parametric.route(port1=(cross0, "P0", 0), port2=(ref_in, "P4", 0)),
        pf.parametric.route(port1=(cross0, "P0", 1), port2=(ref_in, "P2", 1)),
        pf.parametric.route(port1=(cross0, "P0", 2), port2=(ref_in, "P4", 1)),
        pf.parametric.route(port1=(cross0, "P2", 0), port2=(ref_out, "P3", 0)),
        pf.parametric.route(port1=(cross0, "P3", 2), port2=(ref_out, "P5", 1)),
    )

    x_cross += distance + cross_size
    y_cross = -offset / 4
    cross1 = pf.Reference(cross, (x_cross, y_cross), rows=2, spacing=(0, offset / 2))
    c.add(cross1)

    c.add(
        pf.parametric.route_s_bend(port1=(cross0, "P3", i), port2=(cross1, "P1", i))
        for i in range(2)
    )
    c.add(
        pf.parametric.route(port1=(cross0, "P2", i + 1), port2=(cross1, "P0", i))
        for i in range(2)
    )
    c.add(
        pf.parametric.route(port1=(cross1, "P2", 0), port2=(ref_out, "P5", 0)),
        pf.parametric.route(port1=(cross1, "P3", 1), port2=(ref_out, "P3", 1)),
    )

    x_cross += distance + cross_size
    cross2 = pf.Reference(cross, (x_cross, 0))
    c.add(cross2)

    c.add(
        pf.parametric.route_s_bend(port1=(cross1, "P3", 0), port2=(cross2, "P1")),
        pf.parametric.route(port1=(cross1, "P2", 1), port2=(cross2, "P0")),
        pf.parametric.route(port1=(ref_out, "P7", 0), port2=(cross2, "P2")),
        pf.parametric.route(port1=(ref_out, "P1", 1), port2=(cross2, "P3")),
    )

    # Metal pads for the thermal phase-shifters
    num_pads = 34
    pad = pf.Rectangle(size=(pad_size, pad_size))

    # Leave room for routing
    (x_min, y_min), (x_max, y_max) = c.bounds()
    y_min -= 12 * w
    y_max += 9 * w

    # Add all terminals and their structures to the layout
    for i, x in enumerate(
        np.linspace(x_min + pad_size / 2, x_max - pad_size / 2, num_pads)
    ):
        c.add_terminal(
            pf.Terminal("HT", pad.copy().translate((x, y_min - pad_size))),
            f"T{i}",
            add_structure=True,
        )
        c.add_terminal(
            pf.Terminal("HT", pad.copy().translate((x, y_max + pad_size))),
            f"T{i + num_pads}",
            add_structure=True,
        )

    pads = c.terminals

    # Ground net
    c.add(
        pf.parametric.route_manhattan(
            terminal1=pads["T0"],
            terminal2=pads["T34"],
            direction1="y",
            waypoints=[0],
        ),
        pf.parametric.route_manhattan(
            terminal1=(mzi_in, "T0", 0),
            terminal2=pads["T67"],
            direction1="y",
            waypoints=[0],
        ),
        pf.parametric.route_manhattan(
            terminal1=(mzi_in, "T2", 0),
            terminal2=pads["T0"],
            direction1="y",
            waypoints=[0],
        ),
        pf.parametric.route_manhattan(
            terminal1=pads["T33"],
            terminal2=pads["T67"],
            direction1="y",
            waypoints=[0],
        ),
    )

    # Grounded pads
    for ref, t0, t1 in [
        (ref_in, "T0", "T0"),
        (ref_in, "T2", "T2"),
        (ref_in, "T4", "T6"),
        (ref_out, "T1", "T37"),
        (ref_out, "T3", "T39"),
        (ref_out, "T5", "T41"),
        (ref_out, "T7", "T43"),
        (ref_out, "T8", "T44"),
        (ref_out, "T10", "T46"),
        (ref_out, "T48", "T48"),
        (ref_out, "T50", "T50"),
    ]:
        c.add(
            pf.parametric.route_manhattan(
                terminal1=(ref, t0, 0), terminal2=(ref, t1, 1)
            )
        )

    c.add(
        pf.parametric.route_manhattan(
            terminal1=t,
            terminal2=pads[f"T{i + 1}"],
            direction1="y",
            waypoints=[y_min + (20 - i) * w],
        )
        for i, t in enumerate(
            [
                (mzi_in, "T1"),
                (mzi_in, "T3"),
                (ref_in, "T1", 0),
                (ref_in, "T3", 0),
                (ref_in, "T5", 0),
            ]
        )
    )

    # Routes to all phase shifters
    c.add(
        pf.parametric.route_manhattan(
            terminal1=t,
            terminal2=pads[f"T{i + 37}"],
            direction1="y",
            waypoints=[y_max - (18 - i) * w],
        )
        for i, t in enumerate([(ref_in, "T1", 1), (ref_in, "T3", 1), (ref_in, "T7", 1)])
    )

    c.add(
        pf.parametric.route_manhattan(
            terminal1=(ref_in, "T7", 0),
            terminal2=pads["T6"],
            direction1="x",
            waypoints=[ref_in["T7"][0].x_mid + w, y_min + 15 * w],
        ),
        pf.parametric.route_manhattan(
            terminal1=(ref_in, "T5", 1),
            terminal2=pads["T40"],
            direction1="x",
            waypoints=[ref_in["T5"][0].x_mid + w, y_max - 15 * w],
        ),
        pf.parametric.route_manhattan(
            terminal1=(ref_out, "T49", 0), terminal2=pads["T31"], direction1="x"
        ),
        pf.parametric.route_manhattan(
            terminal1=(ref_out, "T51", 0), terminal2=pads["T32"], direction1="y"
        ),
        pf.parametric.route_manhattan(
            terminal1=(ref_out, "T49", 1), terminal2=pads["T65"], direction1="x"
        ),
        pf.parametric.route_manhattan(
            terminal1=(ref_out, "T51", 1), terminal2=pads["T66"], direction1="y"
        ),
    )

    for i in range(4):
        c.add(
            pf.parametric.route_manhattan(
                terminal1=(ref_out, f"T{i * 12}", 0),
                terminal2=pads[f"T{10 - i}"],
                direction1="x",
                waypoints=[ref_out["T0"][0].x_mid - i * w, y_min + (11 + i) * w],
            ),
            pf.parametric.route_manhattan(
                terminal1=(ref_out, f"T{i * 12 + 2}", 0),
                terminal2=pads[f"T{14 - i}"],
                direction1="x",
                waypoints=[ref_out["T2"][0].x_mid - i * w, y_min + (7 + i) * w],
            ),
            pf.parametric.route_manhattan(
                terminal1=(ref_out, f"T{i * 12 + 4}", 0),
                terminal2=pads[f"T{18 - i}"],
                direction1="x",
                waypoints=[ref_out["T4"][0].x_mid - i * w, y_min + (3 + i) * w],
            ),
            pf.parametric.route_manhattan(
                terminal1=(ref_out, f"T{i * 12 + 6}", 0),
                terminal2=pads[f"T{22 - i}"],
                direction1="x",
                waypoints=[ref_out["T6"][0].x_mid - i * w, y_min + max(0, i - 1) * w],
            ),
            pf.parametric.route_manhattan(
                terminal1=(ref_out, f"T{i * 12 + 9}", 0),
                terminal2=pads[f"T{23 + i}"],
                direction1="x",
                waypoints=[
                    ref_out["T9"][0].x_mid + i * w,
                    y_min + abs(i * 2 - 3) // 2 * w,
                ],
            ),
            pf.parametric.route_manhattan(
                terminal1=(ref_out, f"T{i * 12 + 11}", 0),
                terminal2=pads[f"T{27 + i}"],
                direction1="x",
                waypoints=[ref_out["T11"][0].x_mid + i * w, y_min + (1 + i) * w],
            ),
            pf.parametric.route_manhattan(
                terminal1=(ref_out, f"T{i * 12}", 1),
                terminal2=pads[f"T{41 + i}"],
                direction1="x",
                waypoints=[ref_out["T0"][1].x_mid - (3 - i) * w, y_max - (14 - i) * w],
            ),
            pf.parametric.route_manhattan(
                terminal1=(ref_out, f"T{i * 12 + 2}", 1),
                terminal2=pads[f"T{45 + i}"],
                direction1="x",
                waypoints=[ref_out["T2"][1].x_mid - (3 - i) * w, y_max - (10 - i) * w],
            ),
            pf.parametric.route_manhattan(
                terminal1=(ref_out, f"T{i * 12 + 4}", 1),
                terminal2=pads[f"T{49 + i}"],
                direction1="x",
                waypoints=[ref_out["T4"][1].x_mid - (3 - i) * w, y_max - (6 - i) * w],
            ),
            pf.parametric.route_manhattan(
                terminal1=(ref_out, f"T{i * 12 + 6}", 1),
                terminal2=pads[f"T{53 + i}"],
                direction1="x",
                waypoints=[
                    ref_out["T6"][1].x_mid - (3 - i) * w,
                    y_max - (2 - min(i, 2)) * w,
                ],
            ),
            pf.parametric.route_manhattan(
                terminal1=(ref_out, f"T{i * 12 + 9}", 1),
                terminal2=pads[f"T{60 - i}"],
                direction1="x",
                waypoints=[
                    ref_out["T9"][1].x_mid + (3 - i) * w,
                    y_max - abs(3 - 2 * i) // 2 * w,
                ],
            ),
            pf.parametric.route_manhattan(
                terminal1=(ref_out, f"T{i * 12 + 11}", 1),
                terminal2=pads[f"T{64 - i}"],
                direction1="x",
                waypoints=[ref_out["T11"][1].x_mid + (3 - i) * w, y_max - (4 - i) * w],
            ),
        )

    # Make the central ground conductor thicker to minimize heating and inductance
    c.add("HT", pf.Rectangle((pads["T0"].x_mid, 0), (pads["T67"].x_mid, pad_size)))

    c.add_model(pf.CircuitModel())

    return c


device = quantum_processor()
device.write_gds()
viewer(device)