Mach-Zehnder lattice filter

This is an example of a Mach-Zehnder lattice filter with N stages.

import demolib.all as pdk
import ipkiss3.all as i3

import numpy as np
import pylab as plt

The non-dispersive symmetric directional coupler model:

L0=λcarcsin(k0)πΔneffL_0 = \dfrac{\lambda_c arcsin(k_0)}{\pi \Delta n_{eff}}

ϕ0=exp(i2πλnavgLnoncoupling)\phi_0 = exp(i\dfrac{2\pi}{\lambda} n_{avg}L_{non-coupling})

ϕcoupling=exp(i2πλnavgLcoupling)\phi_{coupling} = exp(i\dfrac{2\pi}{\lambda} n_{avg} L_{coupling})

ΔβL=2πλΔneff(Lcoupling+L0)\Delta\beta L = \dfrac{2\pi}{\lambda}\Delta n_{eff} (L_{coupling}+ L_0)

  • Crossing transmission:

Tcross=12(1exp(iΔβL))ϕcouplingϕ0T_{cross} = \dfrac{1}{2}(1 - exp(i\Delta\beta L))\phi_{coupling}\phi_0

  • Straight transmission:

Tstraight=12(1+exp(iΔβL))ϕcouplingϕ0T_{straight} = \dfrac{1}{2}(1 + exp(i\Delta\beta L))\phi_{coupling}\phi_0

class DirCoupModel(i3.CompactModel):
    """
    Non-dispersive symmetric directional coupler model
    """
    parameters = [
        'delta_n',                 # Difference in effective index between the even and odd mode
        'n_avg',                   # Average effective index
        'coupler_length',          # Length of the coupling section
        'non_coupling_length',     # Length of the directional coupler that is not part of the coupling
        'coupling_at_zero_length', # Coupling for a zero length coupling, L_0 = zero length
        'center_wavelength'        # Reference wavelength for all the parameters above
    ]

    terms = [
        i3.OpticalTerm(name='in1'),
        i3.OpticalTerm(name='in2'),
        i3.OpticalTerm(name='out1'),
        i3.OpticalTerm(name='out2'),
    ]

    def calculate_smatrix(parameters, env, S):
        delta_n = parameters.delta_n
        L_0 = parameters.center_wavelength * np.arcsin(parameters.coupling_at_zero_length ) / (np.pi * delta_n)
        phi_0 = np.exp(1j * 2 * np.pi * parameters.n_avg * parameters.non_coupling_length / env.wavelength)
        phi_coupler = np.exp(1j * 2 * np.pi * parameters.n_avg * parameters.coupler_length / env.wavelength)

        delta_beta_x_length = 2 * np.pi / env.wavelength * delta_n * (parameters.coupler_length + L_0)

        # Crossing transmission
        cross = 0.5 * (1.0 - np.exp(1j * delta_beta_x_length)) * phi_coupler * phi_0
        S['in1', 'out2'] = S['out2', 'in1'] = cross
        S['in2', 'out1'] = S['out1', 'in2'] = cross

        # Straight transmission
        straight = 0.5 * (1.0 + np.exp(1j * delta_beta_x_length)) * phi_coupler * phi_0
        S['in1', 'out1'] = S['out1', 'in1'] = straight
        S['in2', 'out2'] = S['out2', 'in2'] = straight

The Directional Coupler described by the previous model:

class DirectionalCoupler(i3.PCell):
    trace_template = i3.WaveguideTemplateProperty(doc="Trace template used for directional coupler")

    def _default_trace_template(self):
        return pdk.SiWireWaveguideTemplate()


    class Layout(i3.LayoutView):
        spacing = i3.PositiveNumberProperty(default=0.2, doc="Spacing between the waveguides", locked=True)
        bend_radius = i3.PositiveNumberProperty(default=10.0, doc="Bend radius of the directional coupler", locked=True)
        coupler_length = i3.PositiveNumberProperty(default=5.0, doc="Length of the straight section of the directional coupler")

        def _generate_instances(self, insts):
            delta = self.trace_template.core_width + self.spacing
            bend_radius = self.bend_radius
            coupler_length = self.coupler_length

            shape = i3.Shape([
                (-coupler_length / 2.0 - bend_radius, bend_radius),
                (-coupler_length / 2.0 - bend_radius, 0.0),
                (-coupler_length / 2.0, 0.0),
                (coupler_length / 2.0, 0.0),
                (coupler_length / 2.0 + bend_radius, 0.0),
                (coupler_length / 2.0 + bend_radius, bend_radius)
            ])

            wg = i3.RoundedWaveguide(name=self.name + "_wg",
                                     trace_template=self.trace_template)
            wg_layout = wg.Layout(shape=shape)

            insts += i3.SRef(reference=wg_layout,
                             name="wav_top",
                             position=(0, delta / 2.0))

            insts += i3.SRef(reference=wg_layout,
                             name="wav_bot",
                             transformation=i3.VMirror() + i3.Translation(translation=(0, -delta / 2.0)))
            return insts

        def _generate_ports(self, ports):
            ports += i3.expose_ports(self.instances, {
                'wav_top:in': 'in2',
                'wav_top:out': 'out2',
                'wav_bot:in': 'in1',
                'wav_bot:out': 'out1',
            })
            return ports


    class CircuitModel(i3.CircuitModelView):
        delta_n_eff = i3.NumberProperty(default=0.04, doc="Difference between even and odd supermode of the dc")
        center_wavelength = i3.PositiveNumberProperty(default=1.31, doc="Reference wavelength for all parameters")
        coupling_at_zero_length = i3.NonNegativeNumberProperty(default=0.03, doc="Field coupling for a zero length coupling")

        def _generate_model(self):
            trace_length = self.layout_view.instances['wav_top'].reference.trace_length()
            coupler_length = self.layout_view.coupler_length
            non_coupling_length = (trace_length - coupler_length)

            return DirCoupModel(delta_n=self.delta_n_eff,
                                n_avg=self.trace_template.n_eff,
                                coupler_length=coupler_length,
                                non_coupling_length=non_coupling_length,
                                coupling_at_zero_length=self.coupling_at_zero_length,
                                center_wavelength=self.center_wavelength
                                )

    class Netlist(i3.NetlistFromLayout):
        pass

The default DirectionalCoupler:

dc = DirectionalCoupler()
dc_layout = dc.Layout()
print("coupler length: {0:.3f}".format(dc_layout.coupler_length))
dc_layout.visualize(annotate=True)
../_images/sphx_glr_plot_mzilattice_001.png

Out:

coupler length: 5.000

The MZI lattice filter that uses an undetermined number (n_stages + 1) of directional couplers.

class MZILattice(i3.PCell):
    """Mach-Zehnder based lattice filter with parameter n_stages.
    """
    directional_couplers = i3.ChildCellListProperty(doc="list of directional couplers", locked=True)
    delay_waveguides = i3.ChildCellListProperty(doc="List of waveguides between the directional couplers")
    n_stages = i3.PositiveIntProperty(doc="Number of stages in the lattice.")
    trace_template = i3.TraceTemplateProperty(doc="The trace template")

    def _default_trace_template(self):
        return pdk.SiWireWaveguideTemplate()

    def _default_directional_couplers(self):
        return [
            DirectionalCoupler(name=self.name + "dc_{}".format(cnt),
                               trace_template=self.trace_template)
            for cnt in range(self.n_stages + 1)
        ]

    def _default_delay_waveguides(self):
        north_wgs = [
            i3.RoundedWaveguide(name=self.name + "wg_north_{}".format(ndx),
                                trace_template=self.trace_template)
            for ndx in range(self.n_stages)]

        south_wgs = [
            i3.RoundedWaveguide(name=self.name + "wg_south_{}".format(ndx),
                                trace_template=self.trace_template)
            for ndx in range(self.n_stages)]

        return north_wgs + south_wgs


    class Layout(i3.LayoutView):
        delay_lengths = i3.ListProperty(default=[100.], doc="List of delay lengths")
        dc_lengths = i3.ListProperty(default=[5., 5.], doc="List of directional coupler lengths")
        bend_radius = i3.PositiveNumberProperty(default=10.0, doc="Radius of the bends")

        def _default_directional_couplers(self):
            dc_models = []
            for idx, dc in enumerate(self.cell.directional_couplers):
                dc_layout = dc.get_default_view(self)
                dc_layout.set(coupler_length=self.dc_lengths[idx])
                dc_models.append(dc_layout)
            return dc_models

        def _default_delay_waveguides(self):
            north_delay_lengths = [max(dl, 0) for dl in self.delay_lengths]
            south_delay_lengths = [abs(min(dl, 0)) for dl in self.delay_lengths]  # negative delay_length -> delay on south side

            min_dl = self.bend_radius
            north_shapes = [i3.Shape([(0.0, 0.0), (0.0, min_dl + dl / 2.0), (2 * self.bend_radius, min_dl + dl / 2.0), (2 * self.bend_radius, 0.0)])
                            for dl in north_delay_lengths
                            ]
            south_shapes = [i3.Shape([(0.0, 0.0), (0.0, -min_dl - dl / 2.0), (2 * self.bend_radius, -min_dl - dl / 2.0), (2 * self.bend_radius, 0.0)])
                            for dl in south_delay_lengths
                            ]

            wgs = []
            for wg, sh in zip(self.cell.delay_waveguides[:self.n_stages], north_shapes):
                wg_lo = wg.get_default_view(i3.LayoutView)
                wg_lo.set(shape=sh, bend_radius=self.bend_radius)
                wgs.append(wg_lo)

            for wg, sh in zip(self.cell.delay_waveguides[self.n_stages:], south_shapes):
                wg_lo = wg.get_default_view(i3.LayoutView)
                wg_lo.set(shape=sh, bend_radius=self.bend_radius)
                wgs.append(wg_lo)

            return wgs

        def _generate_instances(self, insts):
            """Here we join the directional couplers and waveguides with the given delay length and
            place the 'in1' port of the first directional coupler on (0.0, 0.0)."""
            insts0 = {
                "dc_{}".format(cnt): dc
                for cnt, dc in enumerate(self.directional_couplers)
            }

            insts0.update({
                'wg_north_{}'.format(cnt): wg
                for cnt, wg in enumerate(self.delay_waveguides[:self.n_stages])
            })

            insts0.update({
                'wg_south_{}'.format(cnt): wg
                for cnt, wg in enumerate(self.delay_waveguides[self.n_stages:])
            })

            # Determine which ports of the directional couplers and waveguides should be joined:
            south_in_joins = [("dc_{}:out1".format(cnt), "wg_south_{}:in".format(cnt))
                              for cnt in range(self.n_stages)]
            south_out_joins = [("wg_south_{}:out".format(cnt), "dc_{}:in1".format(cnt + 1))
                               for cnt in range(self.n_stages)]
            north_in_joins = [("dc_{}:out2".format(cnt), "wg_north_{}:in".format(cnt))
                              for cnt in range(self.n_stages)]
            north_out_joins = [("wg_north_{}:out".format(cnt), "dc_{}:in2".format(cnt + 1))
                               for cnt in range(self.n_stages)]

            joins = south_in_joins + south_out_joins + north_in_joins + north_out_joins

            insts += i3.place_insts(
                insts=insts0,
                specs=[
                    i3.Place("dc_0:in1", (0.0, 0.0)),
                    i3.Join(joins)
                ]
            )

            return insts

        def _generate_ports(self, ports):
            ports = i3.expose_ports(
                self.instances,
                {"dc_0:in1": "in1",
                 "dc_0:in2": "in2",
                 "dc_{}:out1".format(self.n_stages): "out1",
                 "dc_{}:out2".format(self.n_stages): "out2",
                 }
            )
            return ports

    class Netlist(i3.NetlistFromLayout):
        pass

    class CircuitModel(i3.CircuitModelView):
        delta_n_eff = i3.NumberProperty(default=0.04, doc="Difference between even and odd mode of the dc")
        center_wavelength = i3.PositiveNumberProperty(default=1.31, doc="Reference wavelength for all parameters")
        coupling_at_zero_length = i3.NonNegativeNumberProperty(default=0.03, doc="Field coupling for a zero length coupling")

        def _default_directional_couplers(self):
            dc_models = []
            for dc in self.cell.directional_couplers:
                dc_cm = dc.get_default_view(self)
                dc_cm.set(delta_n_eff=self.delta_n_eff,
                          center_wavelength=self.center_wavelength,
                          coupling_at_zero_length=self.coupling_at_zero_length)
                dc_models.append(dc_cm)
            return dc_models

        def _generate_model(self):
            return i3.HierarchicalModel.from_netlistview(self.netlist_view)

The model parameters:

tt = pdk.SiWireWaveguideTemplate()
tt.Layout().cross_section().visualize()

tt_cm = tt.CircuitModel(n_g=4.18920116149,
                        center_wavelength=1.31,
                        n_eff=2.4579358757)

delta_n_eff = 0.04
coupling_at_zero_length = 0.03
../_images/sphx_glr_plot_mzilattice_002.png

The filter coefficients below lead to a flat-band response.

fsr = 0.04
power_couplings = [0.5, 0.13, 0.12, 0.5, 0.25]
n_eff = tt_cm.n_eff
n_g = tt_cm.n_g
center_wavelength = tt_cm.center_wavelength

L = center_wavelength ** 2 / n_g / fsr
L_pi = center_wavelength / (2 * n_eff)

print("n_eff: {0:.2f}".format(n_eff))
print("n_g:   {0:.2f}".format(n_g))
print("center wavelength: {}".format(center_wavelength))
print("L:    {0:.2f}".format(L))
print("L_pi: {0:.2f}".format(L_pi))

Out:

n_eff: 2.46
n_g:   4.19
center wavelength: 1.31
L:    10.24
L_pi: 0.27

Calculate all layout parameters:

delay_lengths = [L, 2 * L, -(2 * L + L_pi), -2 * L]


def calculate_length_from_coupling(coupling):
    """Calculate the directional coupler length, required to have a certain amplitude coupling."""
    L = center_wavelength * np.arcsin(coupling) / (np.pi * delta_n_eff)
    L0 = center_wavelength * np.arcsin(coupling_at_zero_length) / (np.pi * delta_n_eff)
    return L - L0


dc_lengths = [calculate_length_from_coupling(np.sqrt(kappa)) for kappa in power_couplings]

print("directional coupler lengths = {}".format(dc_lengths))
print("delay lengths = {}".format(delay_lengths))

Out:

directional coupler lengths = [7.8747136068958428, 3.5324806626416718, 3.374845604522557, 7.8747136068958428, 5.1455469402291767]
delay lengths = [10.241212667080566, 20.48242533416113, -20.748909096604063, -20.48242533416113]

Now we can make an MZI lattice with 4 stages, which uses the previously determined trace template. We can visualize it, and write it to a GDSII file.

cell = MZILattice(name="MziLattice", trace_template=tt, n_stages=4)

lv = cell.Layout(delay_lengths=delay_lengths, dc_lengths=dc_lengths)
lv.visualize(annotate=True)
lv.write_gdsii("mzi_lattice.gds")
../_images/sphx_glr_plot_mzilattice_003.png

We can get the S matrix of the circuit model and plot the transmission.

cm = cell.CircuitModel(delta_n_eff=delta_n_eff, center_wavelength=center_wavelength)
wavelengths = np.linspace(1.28, 1.35, 1001)

S = cm.get_smatrix(wavelengths=wavelengths)

plt.figure()
plt.plot(wavelengths, 20 * np.log10(np.abs(S["in1","out1"])), label='Transmission out 1')
plt.plot(wavelengths, 20 * np.log10(np.abs(S["in1","out2"])), label='Transmission out 2')
plt.ylabel("Transmission (dB)")
plt.xlabel("Wavelength ($\mu m$)")
plt.legend()
plt.ylim([-40, 0])
plt.show()
../_images/sphx_glr_plot_mzilattice_004.png

We can now check how sensitive our simulation is with respect to Δneff\Delta n_{eff} . To do this we simulate the same circuit for different values of Δneff\Delta n_{eff} and look at the responses.

plt.figure()

for delta_n_eff in np.linspace(0.03, 0.05, 5):
    cm = cell.CircuitModel(delta_n_eff=delta_n_eff, center_wavelength=center_wavelength)
    wavelengths = np.linspace(1.28, 1.35, 1001)

    S = cm.get_smatrix(wavelengths=wavelengths)

    plt.plot(wavelengths, 20 * np.log10(np.abs(S["in1","out1"])), label='$\Delta n_{eff}$ = '+'{0:.3f}'.format(delta_n_eff)+', $T_{out1}$')
    plt.plot(wavelengths, 20 * np.log10(np.abs(S["in1","out2"])), label='$\Delta n_{eff}$ = '+'{0:.3f}'.format(delta_n_eff)+', $T_{out2}$')
    plt.ylabel("Transmission (dB)")
    plt.xlabel("Wavelength ($\mu m$)")
    plt.legend()
    plt.ylim([-40, 0])

plt.show()
../_images/sphx_glr_plot_mzilattice_005.png