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
from tidy3d.plugins.mode.web import run as run_mode_solver

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("Layers:")
for k, v in sorted(tech.layers.items()):
    print(f"- {k}: {v.layer}")

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

Layers:
- BlackBox: (998, 0)
- Chip design area: (290, 0)
- Deep Trench: (201, 0)
- DevRec: (68, 0)
- Dicing: (210, 0)
- Errors: (999, 0)
- FDTD: (733, 0)
- FbrTgt: (81, 0)
- FloorPlan: (99, 0)
- Keep out: (202, 0)
- M1_heater: (11, 0)
- M2_router: (12, 0)
- M_Open: (13, 0)
- Oxide open (to BOX): (6, 0)
- PinRec: (1, 10)
- PinRecM: (1, 11)
- SEM: (200, 0)
- Si: (1, 0)
- Si N: (20, 0)
- Si N++: (24, 0)
- Si slab: (2, 0)
- SiN: (1, 5)
- Si_Litho193nm: (1, 69)
- Text: (10, 0)
- VC: (40, 0)
- Waveguide: (1, 99)
Ports:
- MM_SiN_TE_1550_3000: 8.0 μm, 6 mode(s)
- MM_TE_1550_2000: 6.0 μm, 10 mode(s)
- MM_TE_1550_3000: 6.0 μm, 14 mode(s)
- Rib_TE_1310_350: 2.35 μm, 1 mode(s)
- Rib_TE_1550_500: 2.5 μm, 1 mode(s)
- SiN_TE_1310_750: 3.0 μm, 1 mode(s)
- SiN_TE_1550_1000: 3.0 μm, 1 mode(s)
- SiN_TE_1550_750: 3.0 μm, 1 mode(s)
- SiN_TE_1550_800: 3.0 μm, 1 mode(s)
- SiN_TE_895_450: 4.0 μm, 2 mode(s)
- SiN_TM_1310_750: 3.0 μm, 2 mode(s)
- SiN_TM_1550_1000: 3.0 μm, 2 mode(s)
- Slot_TE_1550_500: 2.0 μm, 1 mode(s)
- TE-TM_1550_450: 2.0 μm, 2 mode(s)
- TE_1310_350: 2.0 μm, 1 mode(s)
- TE_1310_410: 2.0 μm, 1 mode(s)
- TE_1550_500: 2.0 μm, 1 mode(s)
- TM_1550_500: 2.5 μm, 2 mode(s)
- eskid_TE_1550: 3.31 μm, 1 mode(s)

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]

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(
    port_spec=port_spec,
    coupling_distance=wg_width + coupling_gap,
    coupling_length=1,
    s_bend_offset=offset,
    s_bend_length=6.5 * offset,
    euler_fraction=0.5,
    tidy3d_model_kwargs={"verbose": False},
)

coupler

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"])

Starting...
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()

Starting...
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…
Starting...
Progress: 100%
Testing coupling length = 1.0…
Starting...
Progress: 100%
Testing coupling length = 2.0…
Starting...
Progress: 100%
Testing coupling length = 3.0…
Starting...
Progress: 100%
Testing coupling length = 4.0…
Starting...
Progress: 100%
Testing coupling length = 5.0…
Starting...
Progress: 100%
Testing coupling length = 6.0…
Starting...
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.69, 1.58, 1.37, 5.69, 3.1 ])

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()

Starting...
Progress: 100%
Starting...
Progress: 100%
Starting...
Progress: 100%
Starting...
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 by converting the port specification to a Tidy3D mode solver with to_tidy3d.

mode_solver = tech.ports[port_spec].to_tidy3d(
    [td.C_0 / wavelengths[wavelengths.size // 2]], mesh_refinement=40, group_index=True
)
mode_spec = mode_solver.mode_spec.updated_copy(group_index_step=True)
mode_solver = mode_solver.updated_copy(mode_spec=mode_spec)
mode_data = run_mode_solver(mode_solver)
mode_data.to_dataframe()

08:18:29 EST Mode solver created with
             task_id='fdve-0ef626bc-a81e-4b27-a155-3084630730fb',
             solver_id='mo-f352f239-52dd-46de-8a1b-378a2a15f70c'.

08:18:32 EST Mode solver status: queued

08:19:09 EST Mode solver status: running

08:19:17 EST Mode solver status: success

		wavelength	n eff	k eff	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.409389	0.0	0.982013	0.756527	0.816257	0.200386	4.191854	618.505418

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_data.n_group.isel(f=0, mode_index=0, drop=True))
)
l_pi = wavelengths[wavelengths.size // 2] / (
    2 * float(mode_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.55, 119.11, -118.78, -119.11]

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(
    port_spec=port_spec,
    length=separation,
    name="SEPARATOR",
    waveguide_model_kwargs={"verbose": False},
)

bend = pf.parametric.bend(
    port_spec=port_spec,
    radius=10,
    angle=90,
    euler_fraction=0.5,
    name="BEND",
    waveguide_model_kwargs={"verbose": False},
)


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

    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)

Starting...
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, 0), 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}),
 RandomVariable('bottom_oxide_thickness', **{'value': 3.017, 'stdev': 0.001})]

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': 59.55, 'stdev': 0.029775}),
 RandomVariable('path_difference', **{'value': 118.78, 'stdev': 0.05939}),
 RandomVariable('path_difference', **{'value': 119.11, 'stdev': 0.059555000000000004})]

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,
)

Sample 12 of 12... 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.552638399964366, 119.11893044517554, 118.8577751773862)

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, 0), 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.11", origin=(59.76, 0), rotation=0, scaling=1, x_reflection=False, repetition=Repetition(columns=1, rows=1, spacing=(0, 0))),
 Reference(component="ARMS_119.11", origin=(204.79, 1.37), 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.55', 0),
  {'component_updates': [RandomVariable('path_difference', **{'value': 59.55, 'stdev': 0.029775})]}),
 (('ARMS_119.11', 0),
  {'component_updates': [RandomVariable('path_difference', **{'value': 119.11, 'stdev': 0.059555000000000004})]}),
 (('ARMS_118.78', 0),
  {'component_updates': [RandomVariable('path_difference', **{'value': 118.78, 'stdev': 0.05939})]}),
 (('ARMS_119.11', 1),
  {'component_updates': [RandomVariable('path_difference', **{'value': 119.11, 'stdev': 0.059555000000000004})]})]

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.55', 0, 'DELTA_.*'),
  {'model_updates': [RandomVariable('length', **{'value': 29.775, 'stdev': 0.0148875})]}),
 (('ARMS_119.11', 0, 'DELTA_.*'),
  {'model_updates': [RandomVariable('length', **{'value': 59.555, 'stdev': 0.029777500000000002})]}),
 (('ARMS_118.78', 0, 'DELTA_.*'),
  {'model_updates': [RandomVariable('length', **{'value': 59.39, 'stdev': 0.029695})]}),
 (('ARMS_119.11', 1, 'DELTA_.*'),
  {'model_updates': [RandomVariable('length', **{'value': 59.555, 'stdev': 0.029777500000000002})]})]

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,
)

Sample 12 of 12... 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, 0), xlabel="λ (µm)", ylabel="Transmission (dB)")