optisystem马赫曾德调制器 马赫-曾德

在光子集成线路设计中，马赫曾德干涉仪是一个非常基础和非常重要的器件，其是可编程光子芯片、集成光计算网络以及集成量子密钥分配系统的重要组成部分。接下来介绍利用ipkiss来设计一个可调的MZI，并通过ipkiss来对设计的MZI进行频域和时域的仿真。

老规矩，所有代码如下：

from si_fab import all as pdk
from ipkiss3 import all as i3
import pylab as plt
import numpy as np
import os

class MZI(i3.Circuit):
    """This creates a ppc basic unit, made of tunable MZI.
    Users can specify the length of heater in one of MZI arm，
    They can also specify the length differences between the arms.
    """

    _name_prefix = "MZI"
    dc = i3.ChildCellProperty(doc="Directional coupler in circuit")
    heater = i3.ChildCellProperty(doc="Heater in circuit ")
    waveguide = i3.ChildCellProperty(doc="Waveguide in circuit ")
    wg_buffer_dx = i3.PositiveNumberProperty(default=10.0, doc="Extra non-heated waveguide length around heater")
    heater_length = i3.LockedProperty(doc="Length of heater in arm")
    two_arm_length_difference = i3.NumberProperty(default=0, doc="The length difference between two arm")
    bend_radius = i3.PositiveNumberProperty(default=50, doc="The bend radius of the routing waveguide")
    trace_template = i3.TraceTemplateProperty(
        default=pdk.NWG900(), doc="The trace template for waveguide and heater", locked=True
    )
    v_bias_dc = i3.NumberProperty(default=0.0, doc="DC voltage bias on the MZI (for S-matrix only)")

    def _default_dc(self):
        return pdk.SiNDirectionalCouplerSPower(power_fraction=0.5, target_wavelength=1.55)

    # Set the default value for the heater
    def _default_heater(self):
        ht = pdk.HeatedWaveguide(trace_template=self.trace_template)
        ht.Layout(shape=[(0.0, 0.0), (100.0, 0.0)])
        ht.CircuitModel(v_bias_dc=self.v_bias_dc)
        return ht

    # Set the default value for the parallel waveguide
    def _default_waveguide(self):
        wg = i3.Waveguide(trace_template=self.trace_template)
        wg.Layout(shape=[(0.0, 0.0), (self.heater_length, 0.0)])
        return wg

    def _default_heater_length(self):
        if hasattr(self.heater, "heater_length"):
            return self.heater.heater_length
        else:
            heater_size = self.heater.get_default_view(i3.LayoutView).size_info()
            return heater_size.width

    def _default_insts(self):
        insts = {"dc_1": self.dc, "dc_2": self.dc, "ht": self.heater, "rg": self.waveguide}
        return insts

    def _default_specs(self):
        r = self.bend_radius
        dx = self.wg_buffer_dx
        specs = [
            i3.Place("dc_1:in1", (0.0, 0.0)),
            i3.PlaceRelative("ht:in", "dc_1:out2", (2.0 * r + dx, 2.0 * r)),
            i3.PlaceRelative("rg:in", "dc_1:out1", (2.0 * r + dx, -2.0 * r - self.two_arm_length_difference / 2)),
            i3.PlaceRelative("dc_2:in2", "ht:out", (2.0 * r + dx, -2.0 * r)),
            i3.ConnectManhattan("ht:in", "dc_1:out2", bend_radius=r, min_straight=0),
            i3.ConnectManhattan("rg:in", "dc_1:out1", bend_radius=r, min_straight=0),
            i3.ConnectManhattan("ht:out", "dc_2:in2", bend_radius=r, min_straight=0),
            i3.ConnectManhattan("rg:out", "dc_2:in1", bend_radius=r, min_straight=0),
        ]
        return specs

    def _default_exposed_ports(self):
        exposed_ports = {
            "dc_1:in1": "in1",
            "dc_1:in2": "in2",
            "dc_2:out1": "out1",
            "dc_2:out2": "out2",
            "ht:elec1": "ht1",
            "ht:elec2": "ht2",
        }
        return exposed_ports

def simulate_ppc_unit_linear_voltages(
        cell,
        v_start=0.0,
        v_end=5.0,
        v_num=1000,
        center_wavelength=1.55,
        debug=False,
    ):
    dt = 1e-12  # time step
    t0 = 0e-12  # start time
    t1 = v_num * 1e-12  # end time

    # define optical source
    optical_in = i3.FunctionExcitation(
        port_domain=i3.OpticalDomain,
        excitation_function=lambda t: 1,
    )

    # define electrical source
    def linear_ramp(t):
        return v_end * t / t1 + v_start * (t1 - t) / t1

    electrical_in = i3.FunctionExcitation(
        port_domain=i3.ElectricalDomain,
        excitation_function=linear_ramp,
    )
    gnd = i3.FunctionExcitation(
        port_domain=i3.ElectricalDomain,
        excitation_function=lambda t: 0,
    )
    # Define the child cells and links (direct connections)
    child_cells = {
        "DUT": cell,
        "optical_in": optical_in,
        "ht": electrical_in,
        "gnd": gnd,
        "out1": i3.Probe(port_domain=i3.OpticalDomain),
        "out2": i3.Probe(port_domain=i3.OpticalDomain),
    }
    links = [
        ("optical_in:out", "DUT:in1"),
        ("ht:out", "DUT:ht1"),
        ("gnd:out", "DUT:ht2"),
        ("out1:in", "DUT:out1"),
        ("out2:in", "DUT:out2"),
    ]
    testbench = i3.ConnectComponents(
        child_cells=child_cells,
        links=links,
    )
    # get time simulation results
    cm = testbench.CircuitModel()
    results = cm.get_time_response(
        t0=t0,
        t1=t1,
        dt=dt,
        center_wavelength=center_wavelength,
        debug=debug,
    )
    results.timesteps = results.timesteps[8:]
    results.data = results.data[0:, 8:]
    return results

if __name__ == '__main__':
    # 频域仿真
    mzi = MZI(
        wg_buffer_dx=10.0,
        two_arm_length_difference=0.0,
        bend_radius=50,
        v_bias_dc=0,
    )
    mzi.Layout().visualize(annotate=True)

    # Frequency simulation visualization
    wavelengths = np.linspace(1.5, 1.6, 101)
    cm = mzi.CircuitModel()
    S = cm.get_smatrix(wavelengths=wavelengths)
    plt.plot(wavelengths, i3.signal_power(S["out1", "in1"]), label="out1")
    plt.plot(wavelengths, i3.signal_power(S["out2", "in1"]), label="out2")
    plt.legend()
    plt.xlabel("Wavelengths [um]")
    plt.ylabel("Output optical transmission")
    plt.show()

    # 时域仿真
    # mzi = MZI(
    #     wg_buffer_dx=10.0,
    #     two_arm_length_difference=0.0,
    #     bend_radius=50,
    # )
    # # Time domain simulation
    # # Sweep the voltage for single wavelength
    # results = simulate_ppc_unit_linear_voltages(cell=mzi, v_start=0.0, v_end=5.0, center_wavelength=1.55)
    # times = results.timesteps
    # signals = ["gnd", "ht", "optical_in", "out1", "out2"]
    # ylabels = ["voltage [V]", "voltage [V]", "power [W]", "power [W]", "power [W]"]
    # processes = [np.real, np.real, i3.signal_power, i3.signal_power, i3.signal_power]
    # # Simulation visualization
    # fig, axs = plt.subplots(nrows=len(signals), ncols=1, figsize=(6, 8))
    # for axes, signal, ylabel, process in zip(axs, signals, ylabels, processes):
    #     data = process(results[signal])
    #     axes.set_title(signal)
    #     axes.plot(times, data)
    #     axes.set_ylabel(ylabel)
    #     axes.set_xlabel("time [s]")
    # fig.savefig(os.path.join("{}.png".format("simulation_results_linear_volts")), bbox_inches="tight")
    # plt.show()

代码有点长，先来介绍一下代码的结构：

代码主要分为四个部分：

第一部分：导入ipkiss演示用的pdk包，ipkiss基础包以及python三方包

第二部分：i3.Circuit线路设计函数设计MZI线路

第三部分：定义时域的线路仿真函数

第四部分：MZI实例化以及频域和时域仿真

先来看第二部分，利用i3.Circuit线路设计函数设计MZI线路，这部分代码的结构在教程2中介绍过就不展开了，MZI线路实例化之后如下图：

为了实现MZI的可调，MZI的一个臂加入了热光调制器来调节相位。

在定义MZI属性时，特意加了两个重要的参数：two_arm_length_difference 以及 v_bias_dc

two_arm_length_difference可以方便调节两个臂的长度差，比如：

two_arm_length_difference=200 时

if __name__ == '__main__':
    mzi = MZI(
        wg_buffer_dx=10.0,
        two_arm_length_difference=200.0,
        bend_radius=50,
        v_bias_dc=0,
    )
    mzi.Layout().visualize(annotate=True)

v_bias_dc是调制器两端加的电压，这个参数是为了方便MZI加电压时的频域仿真

第三部分是定义的时域仿真函数

def simulate_ppc_unit_linear_voltages(
        cell,
        v_start=0.0,
        v_end=5.0,
        v_num=1000,
        center_wavelength=1.55,
        debug=False,
    ):

仿真函数的参数主要有：cell是要仿真的线路，电加热器两端的电压v_start，v_end, 仿真电压点数v_num，以及中心波长center_wavelength

仿真函数首先定义仿真时间轴：

dt = 1e-12  # time step
t0 = 0e-12  # start time
t1 = v_num * 1e-12  # end time

定义输入的光源：

# define optical source
    optical_in = i3.FunctionExcitation(
        port_domain=i3.OpticalDomain,
        excitation_function=lambda t: 1,
    )

定义加热器两端的电源：

# define electrical source
    def linear_ramp(t):
        return v_end * t / t1 + v_start * (t1 - t) / t1

    electrical_in = i3.FunctionExcitation(
        port_domain=i3.ElectricalDomain,
        excitation_function=linear_ramp,
    )
    gnd = i3.FunctionExcitation(
        port_domain=i3.ElectricalDomain,
        excitation_function=lambda t: 0,
    )

这里定义的电源正极是随时间线性变化，负极接地（为0）

接下来是定义仿真的单元，这是什么意思呢，相当于模拟一个实际的系统，这个系统需要有哪些器件，比如线路（DUT）,光源（optical_in），电源正极（ht），电源负极（gnd），以及两个探测器（out1，out2），在ipkiss用i3.Probe(port_domain=i3.OpticalDomain)表示虚拟的光探测器

child_cells = {
        "DUT": cell,
        "optical_in": optical_in,
        "ht": electrical_in,
        "gnd": gnd,
        "out1": i3.Probe(port_domain=i3.OpticalDomain),
        "out2": i3.Probe(port_domain=i3.OpticalDomain),
    }

仿真单元定义好，接下来要把它们连起来：

links = [
        ("optical_in:out", "DUT:in1"),
        ("ht:out", "DUT:ht1"),
        ("gnd:out", "DUT:ht2"),
        ("out1:in", "DUT:out1"),
        ("out2:in", "DUT:out2"),
    ]

比如，如果仿真时要在MZI的in1端口输入光，就可以将光源的输出（optical_in:out）连接线路的in1端口（DUT:in1），
而想看线路哪个端口的输出，就可以将探测器的输入（out1:in，out2:in ）接到线路要探测的端口上（DUT:out1，DUT:out2）
电源的连接也是一样(“ht:out”, “DUT:ht1”), (“gnd:out”, “DUT:ht2”)

最后就是是获取时域仿真结果并返回

testbench = i3.ConnectComponents(
        child_cells=child_cells,
        links=links,
    )
    cm = testbench.CircuitModel()
    results = cm.get_time_response(
        t0=t0,
        t1=t1,
        dt=dt,
        center_wavelength=center_wavelength,
        debug=debug,
    )
    results.timesteps = results.timesteps[8:]
    results.data = results.data[0:, 8:]
    return results

第四部分先来看频域仿真

仿真之前，先实例化MZI

mzi = MZI(
        wg_buffer_dx=10.0,
        two_arm_length_difference=0.0,
        bend_radius=50,
        v_bias_dc=0,
    )
    mzi.Layout().visualize(annotate=True)

通过get_smatrix就可以获取整个线路的散射矩阵来进行频域仿真

wavelengths = np.linspace(1.5, 1.6, 101)
cm = mzi.CircuitModel()
S = cm.get_smatrix(wavelengths=wavelengths)
plt.plot(wavelengths, i3.signal_power(S["out1", "in1"]), label="out1")
plt.plot(wavelengths, i3.signal_power(S["out2", "in1"]), label="out2")
plt.legend()
plt.xlabel("Wavelengths [um]")
plt.ylabel("Output optical transmission")
plt.show()

两个端口的输出结果如下：

改变MZI的臂长差（two_arm_length_difference=100.0），以及加电压（v_bias_dc=2）：

mzi = MZI(
        wg_buffer_dx=10.0,
        two_arm_length_difference=100.0,
        bend_radius=50,
        v_bias_dc=2,
    )
mzi.Layout().visualize(annotate=True)
wavelengths = np.linspace(1.5, 1.6, 1001)
cm = mzi.CircuitModel()
S = cm.get_smatrix(wavelengths=wavelengths)
plt.plot(wavelengths, i3.signal_power(S["out1", "in1"]), label="out1")
plt.plot(wavelengths, i3.signal_power(S["out2", "in1"]), label="out2")
plt.legend()
plt.xlabel("Wavelengths [um]")
plt.ylabel("Output optical transmission")
plt.show()

仿真结果如下：

接下来是时域仿真：

仿真之前同样需要先实例化MZI

mzi = MZI(
        wg_buffer_dx=10.0,
        two_arm_length_difference=0.0,
        bend_radius=50,
    )

调用第三部分定义好的时域仿真函数，就能得到MZI时域的响应：

results = simulate_ppc_unit_linear_voltages(cell=mzi, v_start=0.0, v_end=5.0, center_wavelength=1.55)
times = results.timesteps
signals = ["gnd", "ht", "optical_in", "out1", "out2"]
ylabels = ["voltage [V]", "voltage [V]", "power [W]", "power [W]", "power [W]"]
processes = [np.real, np.real, i3.signal_power, i3.signal_power, i3.signal_power]
# Simulation visualization
fig, axs = plt.subplots(nrows=len(signals), ncols=1, figsize=(6, 8))
for axes, signal, ylabel, process in zip(axs, signals, ylabels, processes):
    data = process(results[signal])
    axes.set_title(signal)
    axes.plot(times, data)
    axes.set_ylabel(ylabel)
    axes.set_xlabel("time [s]")
fig.savefig(os.path.join("{}.png".format("simulation_results_linear_volts")), bbox_inches="tight")
plt.show()

仿真结果如下：

可以看到当电源正极的电压随时间线性变化，给线路in1口一个固定的输入，探测器探测到的两个输出端口的光强随时间发生变化。