MZM RF Transmission Line Model (scikit-rf)¶

import copy
from functools import partial

import matplotlib.pyplot as plt
import numpy as np
import skrf as rf
from scipy.constants import speed_of_light as c
from skrf.constants import NumberLike
from skrf.media.definedAEpTandZ0 import DefinedAEpTandZ0

class DefinedAEpTandZ02(DefinedAEpTandZ0):
    """Transmission line model with additional skin-effect losses."""

    def __init__(
        self,
        frequency=None,
        A0=0.0,
        A=0.0,
        f_A=1.0,
        ep_r=1.0,
        tanD=0.0,
        z0_port=None,
        z0=50.0,
        f_low=1.0e6,
        f_high=1.0e12,
        f_ep=1.0e9,
        model="frequencyinvariant",
    ):
        rf.media.Media.__init__(self, frequency=frequency, z0_port=z0_port)
        self.A0 = (A0,)
        self.A, self.f_A = A, f_A
        self.ep_r, self.tanD = ep_r, tanD
        self.z0_characteristic = z0
        self.f_low, self.f_high, self.f_ep = f_low, f_high, f_ep
        self.model = model
        self.z0_port = z0_port

    @property
    def alpha_conductor(self) -> NumberLike:
        """Losses due to conductor resistivity."""
        A, f_A, f, A0 = self.A, self.f_A, self.frequency.f, self.A0
        return A0 + (A * np.log(10) / 20 * np.sqrt(f / f_A))


unloaded_model = partial(
    DefinedAEpTandZ02, f_low=1e9, f_high=100e9, f_ep=20e9, model="djordjevicsvensson"
)

ep_r = 4.04
tanD = 0.11
A = 0.0009
z0 = 85
z0_port = 50
d = 200
units = "um"

freq = rf.Frequency(start=1e-3, stop=40, npoints=100, unit="GHz")
media = unloaded_model(
    frequency=freq, ep_r=ep_r, tanD=tanD, A=A, z0=z0, z0_port=z0_port
)
unloaded_network = media.line(d=d, unit=units)

def vcvs(name, frequency, gain=1):
    """Create a Voltage-Controlled Voltage Source (VCVS) 4-port network."""
    s = rf.ALMOST_ZERO * np.ones([len(frequency), 4, 4], dtype=complex)
    s[:, 0, 0] = 1
    s[:, 1, 1] = 1
    s[:, 2, 0] = gain
    s[:, 2, 1] = -gain
    s[:, 2, 3] = 1
    s[:, 3, 0] = -gain
    s[:, 3, 1] = gain
    s[:, 3, 2] = 1
    return rf.Network(name=name, frequency=frequency, s=s)


def delay_line(frequency, time_delay, zc=50, z0=50, name="delay_line"):
    """Create a transmission line delay 2-port network."""
    gamma = 1j * frequency.w * time_delay
    p = (zc - z0) / (zc + z0)
    prop = np.exp(-gamma)
    denom = 1.0 - (p * prop) ** 2
    s11 = (p * (1.0 - prop**2)) / denom
    s21 = ((1.0 - p**2) * prop) / denom
    s = np.ones([len(frequency), 2, 2], dtype=complex)
    s[:, 0, 0] = s[:, 1, 1] = s11
    s[:, 0, 1] = s[:, 1, 0] = s21
    return rf.Network(name=name, frequency=frequency, s=s)


def voltage_amplifier_2port(frequency, G, Zi, Zo, Z0=50.0, name="amp"):
    """Create a two-port voltage amplifier network."""
    s = np.ones([len(frequency), 2, 2], dtype=complex)
    s[:, 0, 0] = (Zi - Z0) / (Zi + Z0)
    s[:, 0, 1] = 0
    s[:, 1, 0] = 2.0 * G * Zi * Z0 / ((Zi + Z0) * (Zo + Z0))
    s[:, 1, 1] = (Zo - Z0) / (Zo + Z0)
    return rf.Network(name=name, frequency=frequency, s=s)

def generate_loaded_circuit(
    half_unloaded_section,
    r_length=4.7e-3,
    c_per_length=320e-12,
    d=100e-6,
    unit="m",
    z0=50,
    gain=1 / 20.0,
    optical_group_index=4.05,
):
    """Generate a loaded MZM unit cell circuit.

    Args:
        half_unloaded_section: 2-port skrf.Network for half the unloaded line.
        r_length: Sheet resistance of the diode [Ohm*m].
        c_per_length: Capacitance per length of the diode [F/m].
        d: Subsection length.
        unit: Length unit for d.
        z0: Port reference impedance.
        gain: VCVS gain (typically 1/n_sections).
        optical_group_index: Optical group index for the delay line.

    Returns:
        skrf.Circuit: The loaded unit cell circuit.
    """
    half_unloaded_section_in = copy.deepcopy(half_unloaded_section)
    half_unloaded_section_in.name = "unloaded_in"
    half_unloaded_section_out = copy.deepcopy(half_unloaded_section_in)
    half_unloaded_section_out.name = "unloaded_out"

    freq = half_unloaded_section.frequency
    media = DefinedAEpTandZ02(frequency=freq, z0=half_unloaded_section.z0[0, 0])

    d_m = rf.constants.to_meters(d=d, unit=unit)
    cap_per_section = d_m * c_per_length
    r_per_section = r_length / d_m
    time_delay = d_m * optical_group_index / c

    vcvs1 = vcvs(frequency=freq, gain=gain, name="vcvs1")
    delay = delay_line(frequency=freq, time_delay=time_delay, zc=50, name="delay1")
    amp_in = voltage_amplifier_2port(frequency=freq, G=1, Zi=50, Zo=0, name="amp_in")
    amp_out = voltage_amplifier_2port(frequency=freq, G=2, Zi=50, Zo=50, name="amp_out")

    port1 = rf.Circuit.Port(frequency=freq, name="port1", z0=z0)
    port2 = rf.Circuit.Port(frequency=freq, name="port2", z0=z0)
    port3 = rf.Circuit.Port(frequency=freq, name="port3", z0=z0)
    port4 = rf.Circuit.Port(frequency=freq, name="port4", z0=z0)

    cap = media.capacitor(cap_per_section, name="cap1")
    res1 = media.resistor(r_per_section, name="res1")
    ground = rf.Circuit.Ground(frequency=freq, z0=z0, name="ground")

    loaded_connections = [
        [(port1, 0), (half_unloaded_section_in, 0)],
        [(port2, 0), (amp_in, 0)],
        [(port3, 0), (half_unloaded_section_out, 1)],
        [(port4, 0), (amp_out, 1)],
        [(half_unloaded_section_in, 1), (res1, 0), (half_unloaded_section_out, 0)],
        [(res1, 1), (cap, 0), (vcvs1, 1)],
        [(cap, 1), (vcvs1, 0), (ground, 0)],
        [(vcvs1, 2), (amp_in, 1)],
        [(vcvs1, 3), (delay, 0)],
        [(amp_out, 0), (delay, 1)],
    ]

    return rf.Circuit(loaded_connections, name="loaded_section")

n_sections = 10

z0 = 50
section_length = 200.0
section_unit = "um"
freq = rf.Frequency(start=1e-3, stop=40, npoints=100, unit="GHz")
media = unloaded_model(frequency=freq, ep_r=ep_r, tanD=tanD, A=A, z0=60, z0_port=z0)
unloaded_network = media.line(d=0.5 * section_length, unit=section_unit)

unit_cell = generate_loaded_circuit(
    unloaded_network, gain=1.0 / n_sections, d=section_length, unit=section_unit
).network

fig, axes = plt.subplots(ncols=2, figsize=(10, 5))
unit_cell.plot_s_db(3, 0, ax=axes[0])
unit_cell.plot_s_rad(3, 0, ax=axes[1])
plt.suptitle("Loaded unit cell S-parameters")
plt.tight_layout()

n_sections = 10

z0 = 50
zref = 30
section_length = 200.0
section_unit = "um"
freq = rf.Frequency(start=1e-3, stop=40, npoints=100, unit="GHz")
media = unloaded_model(frequency=freq, ep_r=ep_r, tanD=tanD, A=A, z0=60, z0_port=z0)
unloaded_network = media.line(d=0.5 * section_length, unit=section_unit)

unit_cell = generate_loaded_circuit(
    unloaded_network, gain=1.0 / n_sections, d=section_length, unit=section_unit
).network

cascaded_network = rf.network.cascade_list([unit_cell] * n_sections)
amp = voltage_amplifier_2port(frequency=freq, G=1, Zi=50, Zo=0, name="amp")
cascaded_network = rf.connect(cascaded_network, 3, amp, 0)
cascaded_network.name = "Cascaded MZM"
cascaded_network.renormalize([zref, zref, zref, zref])

f_ghz = cascaded_network.frequency.f_scaled
m = cascaded_network.frequency.multiplier
ee_response = cascaded_network.s_db[:, 2, 0]
eo_response = cascaded_network.s_db[:, 3, 0]


def crossing_freq(f, response, threshold):
    """Linear interpolation of the first downward threshold crossing."""
    diff = response - threshold
    idx = np.where(np.diff(np.sign(diff)))[0]
    if len(idx) == 0:
        return np.nan
    i = idx[0]
    t = -diff[i] / (diff[i + 1] - diff[i])
    return f[i] + t * (f[i + 1] - f[i])


f_ee = crossing_freq(f_ghz, ee_response, -6.4)
f_eo = crossing_freq(f_ghz, eo_response, -3.0)

fig, ax = plt.subplots(figsize=(10, 6))
cascaded_network.plot_s_db(2, 0, label="Electrical-Electrical Response")
cascaded_network.frequency.plot(eo_response, label="Electrical-Optical Response")

ax.axhline(y=-6.4, linestyle="--", color="gray", linewidth=0.8)

ax.annotate(
    f"EE BW = {f_ee:.1f} GHz\n($-$6.4 dB)",
    xy=(f_ee * m, -6.4),
    xytext=((f_ee + 1.5) * m, -5.5),
    arrowprops={"arrowstyle": "->", "color": "C0"},
    color="C0",
    fontsize=9,
)
ax.annotate(
    f"EO BW = {f_eo:.1f} GHz\n($-$3 dB)",
    xy=(f_eo * m, -3.0),
    xytext=((f_eo + 1.5) * m, -1.5),
    arrowprops={"arrowstyle": "->", "color": "C1"},
    color="C1",
    fontsize=9,
)

ax.legend()
ax.set_ylim(-9, 1.5)
ax.set_yticks(np.arange(-9, 3.0, 1.5))
ax.set_xlabel("Frequency (GHz)")
ax.set_ylabel("Response (dB)")
ax.set_title(f"MZM EE and EO Response (Zref={zref} Ω)")
plt.tight_layout()