800G Coherent Link Simulation

Download notebook (.zip)

This page simulates the two chips of the series in one link. The link carries dual-polarization 16QAM at 120 GBd. The page ends with the bit error rate (BER) against the optical signal-to-noise ratio (OSNR) for two receiver bandwidths. It compares the BER with the theory curve of the technology page.

The circuit is the layout itself. The two dies of the layout page are instances of one netlist, together with the lasers, 80 km of single-mode fiber and an optical amplifier. Every device of the dies carries the model of its component page:

the FDTD and EME results of the transmitter page and the receiver page, stored as models right after each calculation
an analytic waveguide model on every route, which the layout page sets before it draws the routes

The route lengths, the skew between the X and Y paths and the matched photodiode routes all come from the drawn dies. This page writes the digital signal processing (DSP) of the receiver in numpy, so the reader can follow and change each step.

References

Kikuchi, Kazuro. “Fundamentals of coherent optical fiber communications.” Journal of Lightwave Technology 2016 34 (1), 157-179, doi: 10.1109/JLT.2015.2463719.
Savory, Seb J. “Digital coherent optical receivers: algorithms and subsystems.” IEEE Journal of Selected Topics in Quantum Electronics 2010 16 (5), 1164-1179, doi: 10.1109/JSTQE.2010.2044751.
Pfau, Timo, Sebastian Hoffmann, and Reinhold Noé. “Hardware-efficient coherent digital receiver concept with feedforward carrier recovery for M-QAM constellations.” Journal of Lightwave Technology 2009 27 (8), 989-999, doi: 10.1109/JLT.2008.2010511.

from runpy import run_path

# Coherent_Chip_Layout.qmd defines what this page builds on. Its shared cells are generated into the
# module below; see docs/scripts/generate_notebook_modules.py.
globals().update(
    {k: v for k, v in run_path("misc/coherent_chip_layout.py").items() if not k.startswith("_")}
)

import warnings

from photonforge.live_viewer import LiveViewer
from scipy import signal

# The pads of the receiver die are terminals of the die only, and the circuit leaves them out
warnings.filterwarnings(
    "ignore", message="Terminal .* does not match any reference terminals"
)
# The second input of each hybrid, X, is unused on the die, and no route reaches it
warnings.filterwarnings(
    "ignore", message="These ports are not connected and will be ignored: 'X'"
)

# A viewer for this page; open the printed address in a browser
viewer = LiveViewer()

baud = 120e9  # Symbol rate (Bd)
sps = 4  # Samples per symbol
n_sym = 16384  # Symbols per polarization
dt = 1 / (baud * sps)  # Time step of the simulation (s)
carrier = f0  # Carrier frequency of the time steppers (Hz)
# Frequencies at which each stepper fits its S matrix: 80 % of the band of the time step
fit_freqs = np.linspace(carrier - 0.8 / (2 * dt), carrier + 0.8 / (2 * dt), 81)

levels = np.array([-3, -1, 1, 3]) / np.sqrt(
    10
)  # 16QAM levels per axis, unit average power
gray = np.array([[0, 0], [0, 1], [1, 1], [1, 0]])  # Gray code of each level

LiveViewer started at http://127.0.0.1:35025. Call its 'stop()' method to free the port.

The transmitter die

The transmitter die of the layout page carries everything the link needs:

the models of the transmitter page on its devices, and the model of the MZM on each modulator
a model on every route, and one time stepper for each whole route
the 90 degree bias in the model of each Q heater
the four drive ports XI, XQ, YI and YQ, and a circuit model for the die

# The transmitter die: two fiber facets and four drive ports
print(f"ports of the transmitter die: {sorted(tx_chip.ports)}")
viewer(tx_die)

ports of the transmitter die: ['IN', 'OUT', 'XI', 'XQ', 'YI', 'YQ']

LayoutSchematic

Drive signals

Each of the four tributaries (X and Y, I and Q) carries a random sequence of 16QAM levels. The transmitter DSP shapes them and prepares them for the modulator:

Pulse shaping with a root-raised-cosine (RRC) filter of roll-off 0.1, which confines the signal to $\pm 66$ GHz. The receiver applies the same filter, and the product is a raised cosine with no intersymbol interference.
Digital-to-analog conversion (DAC) with 8 bits and a 60 GHz second-order response, clipped at four standard deviations of the signal.
Pre-distortion by an arcsine. An MZM at its null produces a field proportional to $\sin(\pi V / v_\pi)$ . A drive $V = (v_\pi/\pi)\arcsin(x)$ therefore makes the field proportional to $x$ . The peak is 90 % of full swing.
The electrode response: a filter with the EO response of the 10 mm line from the transmitter page.

These steps are the digital part of the transmitter, so they stay in numpy. The cell then simulates the transmitter die alone, to measure its output power.

v_pi_param = (
    2 * vpil_eff / modulation_length
)  # Voltage that shifts one arm of the MZM by pi (V)
roll = 0.1  # RRC roll-off
dac_bandwidth = 60e9  # Bandwidth of the DAC response (Hz)
dac_bits = 8  # DAC resolution
peak = 0.9  # Peak of the pre-distorted drive, as a fraction of full swing
rng = np.random.default_rng(11)  # Random generator of the transmitted symbols

# Draw the symbol index of each tributary, and the frequency axis of the drive
symbols = {t: rng.integers(0, 4, n_sym) for t in ("XI", "XQ", "YI", "YQ")}
f_axis = np.fft.fftfreq(n_sym * sps, dt)

# Root-raised-cosine magnitude response at the symbol rate
def rrc_spectrum(f):
    af, f1, f2 = np.abs(f), (1 - roll) * baud / 2, (1 + roll) * baud / 2
    # Pass band up to f1, zero above f2
    h = np.where(af <= f1, 1.0, 0.0)
    # Cosine roll-off between f1 and f2
    mid = (af > f1) & (af <= f2)
    h[mid] = np.sqrt(0.5 * (1 + np.cos(np.pi / (roll * baud) * (af[mid] - f1))))
    return h

# DAC response: a second-order Butterworth low-pass, conjugate at negative frequencies
b, a = signal.butter(2, 2 * np.pi * dac_bandwidth, analog=True)
_, h_dac = signal.freqs(b, a, worN=2 * np.pi * np.abs(f_axis))
h_dac = np.where(f_axis >= 0, h_dac, np.conj(h_dac))

# Electrode response of the transmitter page, conjugated to match the sign convention of numpy's FFT
h_eo = np.conj(eo_response(np.abs(f_axis), modulation_length))
h_eo = np.where(f_axis >= 0, h_eo, np.conj(h_eo))

drive = {}
q_dac = 2 ** (dac_bits - 1) - 1  # Largest DAC code
for t, k in symbols.items():
    # Put one impulse per symbol, with the 16QAM level of that symbol
    impulses = np.zeros(n_sym * sps)
    impulses[::sps] = levels[k]
    # Shape the pulses with the RRC filter
    x = np.real(np.fft.ifft(np.fft.fft(impulses) * rrc_spectrum(f_axis)))
    # Quantize to the DAC codes, with full scale at four standard deviations
    full_scale = 4.0 * np.std(x)
    x = np.round(np.clip(x / full_scale, -1, 1) * q_dac) / q_dac
    # Filter with the DAC response
    x = np.real(np.fft.ifft(np.fft.fft(x) * h_dac))
    # Pre-distort by the arcsine, so the field of the MZM follows x
    v = np.arcsin(np.clip(peak * x / np.max(np.abs(x)), -1, 1)) * v_pi_param / np.pi
    # Filter with the electrode response
    drive[t] = np.real(np.fft.ifft(np.fft.fft(v) * h_eo))

# Drive voltages as wave amplitudes on 50 ohm, the input of every simulation below
drive_inputs = pf.TimeSeries(
    {f"{t}@0": (v / np.sqrt(50.0)).astype(complex) for t, v in drive.items()},
    time_step=dt,
)
print(f"peak drive {max(np.max(np.abs(v)) for v in drive.values()):.2f} V")

laser_dbm = 16.0  # Laser power into the input fiber (dBm)
linewidth = 100e3  # Linewidth of the transmitter laser and of the LO (Hz)

# Simulate the transmitter die alone, fed by its laser, to measure its output power
tx_only = pf.component_from_netlist(
    {
        "name": "TX_ONLY",
        "instances": {
            "laser": pfa.cw_laser(
                power=1e-3 * 10 ** (laser_dbm / 10), linewidth=linewidth, seed=3
            ),
            "tx": tx_chip,
        },
        "virtual connections": [(("tx", "IN"), ("laser", "P0"))],
        "ports": [("tx", t, t) for t in ("XI", "XQ", "YI", "YQ")]
        + [("tx", "OUT", "OUT")],
    }
)
tx_only.add_model(pf.CircuitModel(), "Circuit")
with warnings.catch_warnings():
    # The pole fits of the PSR-like devices report errors near 1e-3, which stay small
    warnings.filterwarnings("ignore", message="Fitting error")
    tx_stepper = tx_only.setup_time_stepper(
        dt, carrier, show_progress=False, time_stepper_kwargs={"frequencies": fit_freqs}
    )
tx_out = tx_stepper.step(drive_inputs, show_progress=False)

# Output power in both polarizations, after the start-up of the circuit
startup = 4096  # Samples skipped at the start of each record
e_tx = np.array([np.asarray(tx_out["OUT@0"]), np.asarray(tx_out["OUT@1"])])
tx_out_dbm = 10 * np.log10(
    np.mean(np.sum(np.abs(e_tx[:, startup:]) ** 2, axis=0)) / 1e-3
)
print(f"transmitter output {tx_out_dbm:.2f} dBm in both polarizations")

fig, ax = plt.subplots(1, 1, figsize=(6, 4))
t_ns = np.arange(400) * dt * 1e9
ax.plot(t_ns, drive["XI"][4000:4400], label="XI")
ax.plot(t_ns, drive["XQ"][4000:4400], label="XQ")
ax.set_xlabel("time (ns)")
ax.set_ylabel("drive (V)")
ax.legend()
plt.show()

peak drive 1.75 V
Uploading task 'Mode-LNsinglemoderidgewaveguideforC-bandTEmode…'
Uploading task 'Mode-LNsinglemoderidgewaveguideforC-bandTEmode…'
Uploading task 'Mode-LNstripwaveguideforC-bandTEmode…'
Uploading task 'Mode-LNstripwaveguideforC-bandTEmode…'
Uploading task 'Mode-LNsinglemoderidgewaveguideforC-bandTEmode…'
Uploading task 'Mode-LNsinglemoderidgewaveguideforC-bandTEmode…'
Starting task 'Mode-LNsinglemoderidgewaveguideforC-bandTEmode': https://tidy3d.simulation.cloud/workbench?taskId=mo-04474ed2-b771-4f12-ba2c-a2f67ffc9f98
Starting task 'Mode-LNsinglemoderidgewaveguideforC-bandTEmode': https://tidy3d.simulation.cloud/workbench?taskId=mo-feb3805b-53fa-4f57-9d5c-24f6bbdc9a9d
Starting task 'Mode-LNstripwaveguideforC-bandTEmode': https://tidy3d.simulation.cloud/workbench?taskId=mo-cf5fffae-183d-4392-848b-f364bc891740
Starting task 'Mode-LNsinglemoderidgewaveguideforC-bandTEmode': https://tidy3d.simulation.cloud/workbench?taskId=mo-949c727b-065f-407d-b575-a639a7f25df2
Starting task 'Mode-LNstripwaveguideforC-bandTEmode': https://tidy3d.simulation.cloud/workbench?taskId=mo-4a5560fe-c1b9-49fc-95db-248cf492041c
Downloading data from 'Mode-LNsinglemoderidgewaveguideforC-bandTEmode'…
Starting task 'Mode-LNsinglemoderidgewaveguideforC-bandTEmode': https://tidy3d.simulation.cloud/workbench?taskId=mo-88bdedf8-9524-42dd-84f9-7312124ceeb3
Downloading data from 'Mode-LNsinglemoderidgewaveguideforC-bandTEmode'…
Downloading data from 'Mode-LNstripwaveguideforC-bandTEmode'…
Downloading data from 'Mode-LNsinglemoderidgewaveguideforC-bandTEmode'…
Downloading data from 'Mode-LNstripwaveguideforC-bandTEmode'…
Downloading data from 'Mode-LNsinglemoderidgewaveguideforC-bandTEmode'…
transmitter output -1.00 dBm in both polarizations

The transmitter delivers about $-1$ dBm from a 16 dBm laser. The modulators sit at their null, and the 16QAM field averages well below the peak the modulator can deliver. The edge couplers, the MMIs, the modulator arms and the routes add several decibels more.

The receiver die

The receiver die of the layout page carries its models too. Each photodiode is its own cell, with the front-end model of the receiver page and its own noise seed. Its output E0 is a port of the die, named after the hybrid output that feeds it. coherent_rx_die takes the two bandwidths of the front end and a seed. Each simulation below therefore builds the die again, with new noise.

The page compares two front ends:

A: the 67 GHz photodiode with the 80 GHz TIA of the receiver page, 47 GHz together
B: two equal poles at 124 GHz, 80 GHz together, the bandwidth of the 128 GBd receivers on the market

fronts = {
    "A": (67e9, 80e9),
    "B": (124.3e9, 124.3e9),
}  # Photodiode and TIA bandwidths (Hz)
outputs = ("Ip", "In", "Qp", "Qn")  # Hybrid outputs of each polarization
pd_names = [
    pol + o for pol in ("X", "Y") for o in outputs
]  # Front-end outputs of the receiver die

# The receiver die: the two fiber inputs and the eight front-end outputs
print(f"ports of the receiver die: {sorted(rx_die.ports)}")
viewer(rx_die)

ports of the receiver die: ['LO', 'SIG', 'XIn', 'XIp', 'XQn', 'XQp', 'YIn', 'YIp', 'YQn', 'YQp']

The link circuit

The cell below builds one netlist for each simulation. Both cases start with the transmitter laser, the transmitter die and the receiver die. The LO is a second laser with the same 100 kHz linewidth, 1.5 GHz above the transmitter laser, as two free-running lasers are. The DSP must find and remove that offset.

Back to back: a variable optical attenuator (VOA) between the dies sets $-5$ dBm at the receiver fiber.
80 km link: a VOA sets the launch power, followed by the fiber, its polarization mode dispersion (PMD) and an erbium-doped fiber amplifier (EDFA).

The parts of the link come from the abstract library:

Fiber: fiber with both polarizations, 80 km with 0.2 dB/km of loss and its chromatic dispersion.
PMD: pmd_fiber, a random realization of the fiber’s birefringence with a mean differential group delay of 1 ps.
EDFA: one optical_amplifier per polarization, between two polarization_splitter_rotator blocks, with a 5 dB noise figure. Its gain control holds the receiver input at $-5$ dBm.
VOA: 1 m of fiber with the loss of the attenuator and no dispersion.
LO input: a polarization splitter-rotator puts the LO into the TE mode of the receiver fiber. A single-mode port connected to a two-mode port would pass no light. The unused input of the splitter-rotator ends in a termination.

The launch power sets the OSNR. The amplifier gain $G$ is the 16 dB span loss plus the step from the launch power to $-5$ dBm. The ASE of both polarizations in 0.1 nm is

P_\text{ASE} = 2\, n_\text{sp} (G - 1)\, h\nu\, B_\text{ref}, \qquad n_\text{sp} = \frac{F}{2},

with $F$ the noise figure. The OSNR is the $-5$ dBm of signal over $P_\text{ASE}$ .

p_rx_dbm = -5.0  # Signal power at the receiver fiber (dBm)
lo_dbm = 0.0  # LO power at its fiber (dBm)
lo_offset = 1.5e9  # LO frequency above the transmitter laser (Hz)
length_km, alpha_db_km = 80.0, 0.2  # Fiber length (km) and loss (dB/km)
span_db = length_km * alpha_db_km  # Fiber loss (dB)
nf_db = 5.0  # EDFA noise figure (dB)
mean_dgd = 1.0  # Mean differential group delay of the fiber (ps)
voa_km = 1e-3  # Length of the fiber that stands for a VOA (km)
h_nu = 6.62607015e-34 * carrier  # Photon energy (J)

# The 80 km fiber, with both polarizations; the DSP also uses it to compensate the dispersion
fiber = pfa.fiber(length_km=length_km, attenuation=alpha_db_km, dual_polarization=True)

# VOA: a short fiber with the loss of the attenuator and no dispersion
def voa(loss_db):
    return pfa.fiber(
        length_km=voa_km,
        attenuation=loss_db / voa_km,
        dispersion=0,
        dispersion_slope=0,
        dual_polarization=True,
    )

# The link circuit: back to back if launch_dbm is None, otherwise the 80 km link at that launch power
def build_link(front, seed, launch_dbm=None):
    f_pd, f_tia = fronts[front]
    # Lasers, dies and the splitter-rotator that puts the LO into the TE mode (mode 1) of its fiber
    instances = {
        "laser": pfa.cw_laser(
            power=1e-3 * 10 ** (laser_dbm / 10), linewidth=linewidth, seed=3
        ),
        "tx": tx_chip,
        "rx": coherent_rx_die(f_pd=f_pd, f_tia=f_tia, seed=seed),
        "lo": pfa.cw_laser(
            power=1e-3 * 10 ** (lo_dbm / 10),
            linewidth=linewidth,
            seed=seed + 1000,
            frequency=carrier + lo_offset,
        ),
        "lo_pol": pfa.polarization_splitter_rotator(),
        "lo_end": pfa.optical_termination(),
    }
    # The LO enters the splitter-rotator at P2, and its unused P1 ends in a termination
    connections = [
        (("tx", "IN"), ("laser", "P0")),
        (("lo_pol", "P2"), ("lo", "P0")),
        (("lo_end", "P0"), ("lo_pol", "P1")),
        (("rx", "LO"), ("lo_pol", "P0")),
    ]
    if launch_dbm is None:
        # Back to back: one VOA between the dies
        instances["voa"] = voa(tx_out_dbm - p_rx_dbm)
        connections += [(("voa", "P0"), ("tx", "OUT")), (("rx", "SIG"), ("voa", "P1"))]
    else:
        # EDFA gain: the span loss plus the step from the launch power to the receiver input
        gain_db = span_db + p_rx_dbm - launch_dbm
        instances |= {
            "voa": voa(tx_out_dbm - launch_dbm),
            "fiber": fiber,
            "pmd": pfa.pmd_fiber(
                length_km=length_km, pmd=mean_dgd / np.sqrt(length_km), seed=100 + seed
            ),
            "split": pfa.polarization_splitter_rotator(),
            "edfa_0": pfa.optical_amplifier(
                gain=gain_db, noise_figure=nf_db, seed=200 + 10 * seed
            ),
            "edfa_1": pfa.optical_amplifier(
                gain=gain_db, noise_figure=nf_db, seed=201 + 10 * seed
            ),
            "combine": pfa.polarization_splitter_rotator(),
        }
        # VOA, fiber, PMD, one EDFA per polarization, and the receiver
        connections += [
            (("voa", "P0"), ("tx", "OUT")),
            (("fiber", "P0"), ("voa", "P1")),
            (("pmd", "P0"), ("fiber", "P1")),
            (("split", "P0"), ("pmd", "P1")),
            (("edfa_0", "P0"), ("split", "P1")),
            (("edfa_1", "P0"), ("split", "P2")),
            (("combine", "P1"), ("edfa_0", "P1")),
            (("combine", "P2"), ("edfa_1", "P1")),
            (("rx", "SIG"), ("combine", "P0")),
        ]
    # Ports: the four drive inputs and the eight photodiode outputs
    ports = [("tx", t, t) for t in ("XI", "XQ", "YI", "YQ")] + [
        ("rx", n, n) for n in pd_names
    ]
    link = pf.component_from_netlist(
        {
            "name": "LINK",
            "instances": instances,
            "virtual connections": connections,
            "ports": ports,
        }
    )
    link.add_model(pf.CircuitModel(), "Circuit")
    return link

# OSNR in 0.1 nm at the EDFA output for a launch power
def osnr_db(launch_dbm):
    gain = 10 ** ((span_db + p_rx_dbm - launch_dbm) / 10)
    p_ase = 2 * 10 ** (nf_db / 10) / 2 * (gain - 1) * h_nu * b_ref
    return 10 * np.log10(1e-3 * 10 ** (p_rx_dbm / 10) / p_ase)

# Balanced outputs in volts: I = Ip - In, and Q = Qn - Qp, which gives the field rather than its conjugate
def balanced(out, pol):
    v = {o: np.real(np.asarray(out[f"{pol}{o}@0"])) * np.sqrt(50.0) for o in outputs}
    return v["Ip"] - v["In"], v["Qn"] - v["Qp"]

# 8-bit analog-to-digital converter, full scale at four standard deviations
def adc(x, bits=8):
    full = 4.0 * np.std(x)
    q = 2 ** (bits - 1) - 1
    return np.round(np.clip(x / full, -1, 1) * q) / q * full

Receiver DSP

The four ADC outputs go through the standard chain of a coherent receiver [2]:

IQ orthogonalization. The Gram-Schmidt procedure removes the part of Q that correlates with I. This corrects the phase error of the hybrid. The procedure then normalizes both.
Chromatic dispersion compensation. The dispersion of the fiber is a known all-pass phase. The DSP applies the inverse phase in the frequency domain, together with the RRC matched filter, and keeps two samples per symbol.
Frequency offset estimation. The fourth power of 16QAM removes most of its modulation and leaves a tone at four times the offset between the lasers.
Polarization demultiplexing. A 2x2 multiple-input multiple-output (MIMO) equalizer with 21 taps at half-symbol spacing undoes the polarization rotation and the PMD of the fiber. It also equalizes the bandwidth of the chips. The least-mean-squares (LMS) algorithm trains its taps on a known preamble of 4096 symbols, as a real frame carries. The taps then stay constant, because the channel does not change over the 136 ns record.
Carrier phase recovery. Blind phase search [3] tests 32 phases per quarter turn. It keeps the phase whose decisions sit closest to the constellation, averaged over 128 symbols. This tracks the phase noise of the two lasers.
Pilots. Blind phase search leaves a quarter-turn ambiguity. One known pilot symbol in every 64 sets the quadrant, with a vote over five neighbouring pilots.
Decisions and BER. Gray decoding against the transmitted bits, after 3000 symbols for the equalizer to settle.

def dsp(
    outs,
    fiber=None,
    n_taps=21,
    bps_window=128,
    n_train=4096,
    train_mus=(2e-3, 5e-4, 1e-4),
):
    # 1) and 2): orthogonalize I and Q, compensate the dispersion, apply the matched filter
    f4 = np.fft.fftfreq(outs.shape[1], dt)
    z = []
    for p in range(2):
        # Remove the mean of I and Q, and normalize I
        i, q = (
            outs[2 * p] - outs[2 * p].mean(),
            outs[2 * p + 1] - outs[2 * p + 1].mean(),
        )
        i = i / np.sqrt(np.mean(i**2))
        # Remove the part of Q that correlates with I, and normalize Q
        q = q - np.mean(q * i) * i
        q = q / np.sqrt(np.mean(q**2))
        # Inverse phase of the fiber; back to back there is nothing to undo
        cd = (
            1.0
            if fiber is None
            else np.exp(-1j * np.angle(fiber.s_matrix(carrier - f4)[("P0@0", "P1@0")]))
        )
        # Apply the inverse phase and the RRC filter, then keep two samples per symbol
        zz = np.fft.ifft(np.fft.fft(i + 1j * q) * cd * rrc_spectrum(f4))
        z.append(zz[:: sps // 2])
    z = np.array(z)
    # Normalize each polarization to unit power
    z /= np.sqrt(np.mean(np.abs(z) ** 2, axis=1, keepdims=True))

    # 3) Frequency offset: the peak of the 4th-power spectrum sits at 4 times the offset
    n = z.shape[1]
    spectrum = np.abs(np.fft.fft(z[0] ** 4, 8 * n))
    f_offset = np.fft.fftfreq(8 * n, 1 / (2 * baud))[int(np.argmax(spectrum))] / 4
    # Remove the offset from both polarizations
    z *= np.exp(-2j * np.pi * f_offset * np.arange(n) / (2 * baud))

    # 4) Frame sync: find the sample phase and the lag where the preamble starts
    ref = {
        pol: levels[symbols[pol + "I"]] + 1j * levels[symbols[pol + "Q"]]
        for pol in ("X", "Y")
    }
    best = None
    for phase in (0, 1):
        for p in range(2):
            for pol in ("X", "Y"):
                # Cross-correlate one sample phase of one output with one transmitted polarization
                zz = z[p, phase::2]
                m = min(len(zz), n_sym)
                cc = np.abs(
                    np.fft.ifft(
                        np.fft.fft(zz[:m], 2 * m)
                        * np.conj(np.fft.fft(ref[pol][:m], 2 * m))
                    )
                )[:m]
                if best is None or cc.max() > best[0]:
                    best = (cc.max(), int(np.argmax(cc)), phase)
    _, lag0, phase0 = best
    # Start the equalizer as a pass-through: a unit centre tap on each diagonal
    w = np.zeros((2, 2, n_taps), complex)
    w[0, 0, n_taps // 2] = w[1, 1, n_taps // 2] = 1
    ns = n // 2 - n_taps  # Number of equalizer outputs
    # Output k carries the transmitted symbol k - shift
    shift = lag0 - (n_taps // 2 - phase0) // 2
    k_train = np.arange(max(0, shift), min(ns, shift + n_train + max(0, -shift)))
    # Train the taps with LMS on the preamble, with a decreasing step
    for mu in train_mus:
        for k in k_train:
            u = z[:, 2 * k : 2 * k + n_taps][:, ::-1]
            yk = np.einsum("abt,bt->a", w, u)
            d = np.array([ref["X"][k - shift], ref["Y"][k - shift]])
            w -= mu * (yk - d)[:, None, None] * np.conj(u)[None, :, :]
    # Equalize the whole record with the fixed taps
    y = np.zeros((2, ns), complex)
    for k in range(ns):
        y[:, k] = np.einsum("abt,bt->a", w, z[:, 2 * k : 2 * k + n_taps][:, ::-1])

    # 5) Blind phase search, and 6) the quadrant from the pilots
    const = (
        levels[:, None] + 1j * levels[None, :]
    ).ravel()  # The 16 constellation points
    phis = np.arange(32) / 32 * np.pi / 2  # Test phases over a quarter turn
    out = np.zeros_like(y)
    kk = np.arange(ns)
    for p, pol in enumerate(("X", "Y")):
        # Distance of each rotated symbol to its nearest constellation point, averaged over the window
        rot = y[p][None, :] * np.exp(1j * phis)[:, None]
        dist = np.min(np.abs(rot[:, :, None] - const[None, None, :]) ** 2, axis=2)
        dist = signal.fftconvolve(dist, np.ones((1, bps_window)), mode="same", axes=1)
        # Pick the best test phase per symbol, unwrapped with a quarter-turn period
        phase = np.unwrap(phis[np.argmin(dist, axis=0)] * 4) / 4
        yc = y[p] * np.exp(1j * phase)
        # Pilots: one known symbol in every 64 of the transmitted sequence
        m_ref = kk - shift
        pilots = kk[(m_ref >= 0) & (m_ref < n_sym) & (m_ref % 64 == 0)]
        # Quadrant of each pilot, then a majority vote over five neighbouring pilots
        q = np.mod(
            np.round(
                np.angle(ref[pol][pilots - shift] * np.conj(yc[pilots])) / (np.pi / 2)
            ),
            4,
        )
        votes = np.array(
            [
                np.bincount(q[max(0, i - 2) : i + 3].astype(int), minlength=4).argmax()
                for i in range(len(q))
            ]
        )
        # Apply the quadrant of the nearest pilot to each symbol
        quadrant = votes[np.clip(np.searchsorted(pilots, kk - 32), 0, len(pilots) - 1)]
        out[p] = yc * 1j ** np.round(quadrant)

    # 7) Assign each output to its polarization, decide, and count bit errors
    res = {}
    # Skip the settling of the equalizer at the start, and the wrap of the circular filter at the end
    skip, tail = 3000, 1000
    for p in range(2):
        # Find the transmitted polarization, the quarter turn and the lag that match this output best
        best = None
        for pol in ("X", "Y"):
            for r4 in range(4):
                o = out[p] * 1j**r4
                cc = np.abs(
                    np.fft.ifft(
                        np.fft.fft(o, 2 * ns) * np.conj(np.fft.fft(ref[pol], 2 * ns))
                    )
                )[:ns]
                lag = int(np.argmax(cc))
                if best is None or cc[lag] > best[0]:
                    best = (cc[lag], pol, r4, lag)
        _, pol, r4, lag = best
        # Align the output with the transmitted symbols and remove the complex gain
        o = (out[p] * 1j**r4)[lag + skip : ns - tail]
        r = ref[pol][skip : ns - lag - tail]
        o = o / (np.vdot(r, o) / np.vdot(r, r))
        # Error vector magnitude against the transmitted symbols
        evm = np.sqrt(np.mean(np.abs(o - r) ** 2) / np.mean(np.abs(r) ** 2))
        # Decide each axis to the nearest level and count the Gray-coded bit errors
        decide = lambda x: np.argmin(np.abs(x[:, None] - levels[None, :]), axis=1)
        errors = np.sum(gray[decide(o.real)] != gray[decide(r.real)]) + np.sum(
            gray[decide(o.imag)] != gray[decide(r.imag)]
        )
        res[pol] = dict(
            evm=float(evm), errors=int(errors), bits=int(4 * len(r)), symbols=o
        )
    return res, float(f_offset)

# Simulate one link circuit, detect, and run the DSP
def simulate(front, seed, launch_dbm=None):
    link = build_link(front, seed, launch_dbm)
    with warnings.catch_warnings():
        # The pole fits of the PSR-like devices report errors near 1e-3, which stay small
        warnings.filterwarnings("ignore", message="Fitting error")
        stepper = link.setup_time_stepper(
            dt,
            carrier,
            show_progress=False,
            time_stepper_kwargs={"frequencies": fit_freqs},
        )
    out = stepper.step(drive_inputs, show_progress=False)
    # Digitize the four balanced outputs: XI, XQ, YI and YQ
    outs = np.array([adc(v) for pol in ("X", "Y") for v in balanced(out, pol)])
    # Run the DSP, with the dispersion compensation only on the 80 km link
    res, f_offset = dsp(outs, None if launch_dbm is None else fiber)
    # BER over both polarizations
    ber = sum(r["errors"] for r in res.values()) / sum(r["bits"] for r in res.values())
    return ber, res, f_offset

Back to back

The first simulation connects the transmitter die straight to the receiver die through the VOA. What remains are the impairments of the two chips, the lasers and the converters:

the limited bandwidth of the modulators and the front end
the phase error and the imbalance of the hybrids
the polarization-dependent losses of the couplers and the rotators
the quantization of the DAC and the ADC
the noise of the photodiodes and the TIAs

# Simulate back to back with front end A
ber_b2b, res, f_offset = simulate("A", seed=0)
print(f"estimated laser offset {f_offset * 1e-9:.3f} GHz")
for pol in ("X", "Y"):
    print(
        f"{pol}: EVM {100 * res[pol]['evm']:.2f} %, {res[pol]['errors']} errors in {res[pol]['bits']} bits"
    )

fig, ax = plt.subplots(1, 2, figsize=(9, 4.5))
for a, pol in zip(ax, ("X", "Y")):
    o = res[pol]["symbols"]
    a.plot(o.real, o.imag, ".", ms=1, alpha=0.3)
    a.set_aspect("equal")
    a.set_xlim(-1.5, 1.5)
    a.set_ylim(-1.5, 1.5)
    a.set_title(f"{pol} polarization, back to back")
plt.show()

Uploading task 'Mode-StripwaveguideTE500nm…'
Uploading task 'Mode-StripwaveguideTE500nm…'
Uploading task 'Mode-StripwaveguideTE500nm…'
Uploading task 'Mode-StripwaveguideTE500nm…'
Uploading task 'Mode-StripwaveguideTE500nm…'
Uploading task 'Mode-StripwaveguideTE500nm…'
Starting task 'Mode-StripwaveguideTE500nm': https://tidy3d.simulation.cloud/workbench?taskId=mo-5d42e4f0-c24c-495d-ab25-70127c5a1ef2
Starting task 'Mode-StripwaveguideTE500nm': https://tidy3d.simulation.cloud/workbench?taskId=mo-8fdc9520-f395-4268-b500-26b13395b80f
Starting task 'Mode-StripwaveguideTE500nm': https://tidy3d.simulation.cloud/workbench?taskId=mo-7509df48-95a4-415c-953e-8bc5e50d5cca
Starting task 'Mode-StripwaveguideTE500nm': https://tidy3d.simulation.cloud/workbench?taskId=mo-6cd007c1-7e2a-4e67-98eb-3de6a274686c
Starting task 'Mode-StripwaveguideTE500nm': https://tidy3d.simulation.cloud/workbench?taskId=mo-80f61838-4b77-4905-ac05-a8cf5fe14a24
Downloading data from 'Mode-StripwaveguideTE500nm'…
Starting task 'Mode-StripwaveguideTE500nm': https://tidy3d.simulation.cloud/workbench?taskId=mo-348fa8da-1d4a-4455-a855-203e9c07704f
Downloading data from 'Mode-StripwaveguideTE500nm'…
Downloading data from 'Mode-StripwaveguideTE500nm'…
Downloading data from 'Mode-StripwaveguideTE500nm'…
Downloading data from 'Mode-StripwaveguideTE500nm'…
Downloading data from 'Mode-StripwaveguideTE500nm'…
estimated laser offset 1.499 GHz
X: EVM 8.84 %, 0 errors in 49396 bits
Y: EVM 8.86 %, 0 errors in 49396 bits

The back-to-back constellations are clean, with an error vector magnitude (EVM) of about 9 % and no bit errors. The noise of the TIAs and the quantization of the converters set most of the EVM.

The 80 km link

The cell below sweeps the launch power from $-18$ to $-2$ dBm for both front ends. The OSNR then goes from about 19 to 35 dB. Each point is a new circuit with a new realization of the PMD. It also has new seeds for the noise of the lasers, the amplifiers and the photodiodes.

launches = [
    -18.0,
    -16.0,
    -14.0,
    -12.0,
    -10.0,
    -8.0,
    -4.0,
    -2.0,
]  # Launch powers into the fiber (dBm)
results = {}
for front in fronts:
    rows = []
    for k, launch in enumerate(launches):
        # Simulate the link at this launch power, with seed k for every random process
        ber, res, _ = simulate(front, seed=k, launch_dbm=launch)
        rows.append((osnr_db(launch), ber, res))
        print(
            f"front end {front}, launch {launch:+5.1f} dBm: OSNR {osnr_db(launch):5.1f} dB, BER {ber:.2e}, "
            f"EVM X {100 * res['X']['evm']:.1f} % Y {100 * res['Y']['evm']:.1f} %"
        )
    results[front] = rows

front end A, launch -18.0 dBm: OSNR  19.0 dB, BER 1.06e-01, EVM X 39.4 % Y 45.4 %
front end A, launch -16.0 dBm: OSNR  21.0 dB, BER 6.04e-02, EVM X 30.0 % Y 33.8 %
front end A, launch -14.0 dBm: OSNR  23.0 dB, BER 2.92e-02, EVM X 24.1 % Y 26.5 %
front end A, launch -12.0 dBm: OSNR  25.0 dB, BER 1.29e-02, EVM X 20.0 % Y 21.8 %
front end A, launch -10.0 dBm: OSNR  27.0 dB, BER 3.74e-03, EVM X 16.3 % Y 18.4 %
front end A, launch  -8.0 dBm: OSNR  29.0 dB, BER 9.50e-04, EVM X 14.1 % Y 15.3 %
front end A, launch  -4.0 dBm: OSNR  33.1 dB, BER 5.52e-05, EVM X 11.8 % Y 12.3 %
front end A, launch  -2.0 dBm: OSNR  35.2 dB, BER 1.10e-05, EVM X 11.5 % Y 10.9 %
front end B, launch -18.0 dBm: OSNR  19.0 dB, BER 1.02e-01, EVM X 38.5 % Y 43.1 %
front end B, launch -16.0 dBm: OSNR  21.0 dB, BER 5.78e-02, EVM X 29.5 % Y 33.2 %
front end B, launch -14.0 dBm: OSNR  23.0 dB, BER 2.68e-02, EVM X 23.6 % Y 26.0 %
front end B, launch -12.0 dBm: OSNR  25.0 dB, BER 1.07e-02, EVM X 19.4 % Y 21.1 %
front end B, launch -10.0 dBm: OSNR  27.0 dB, BER 2.68e-03, EVM X 15.8 % Y 17.4 %
front end B, launch  -8.0 dBm: OSNR  29.0 dB, BER 4.53e-04, EVM X 13.2 % Y 14.4 %
front end B, launch  -4.0 dBm: OSNR  33.1 dB, BER 0.00e+00, EVM X 10.3 % Y 11.1 %
front end B, launch  -2.0 dBm: OSNR  35.2 dB, BER 0.00e+00, EVM X 9.5 % Y 10.1 %

BER against OSNR

The cell below plots both front ends against the theory curve of the technology page. It also finds the OSNR that each front end needs at the FEC threshold, $2 \times 10^{-2}$ , and at $10^{-3}$ . It interpolates $\log_{10}$ BER between the simulated points.

osnr_axis = np.linspace(18, 38, 200)  # OSNR axis of the theory curve (dB)
ber_theory = ber_16qam(10 ** (osnr_axis / 10) * b_ref / baud)

fig, ax = plt.subplots(1, 1, figsize=(6.5, 4.5))
ax.semilogy(osnr_axis, ber_theory, "k-", label="theory")
labels = {"A": "47 GHz front end", "B": "80 GHz front end"}
for front, rows in results.items():
    osnr = np.array([r[0] for r in rows])
    ber = np.array([max(r[1], 1e-6) for r in rows])
    ax.semilogy(osnr, ber, "o-", label=labels[front])
    # Interpolate log BER on the points with errors, for each target BER
    nonzero = np.array([r[1] for r in rows]) > 0
    for target in (ber_fec, 1e-3):
        need = np.interp(-np.log10(target), -np.log10(ber[nonzero]), osnr[nonzero])
        print(f"{labels[front]}: required OSNR {need:.1f} dB at BER {target:.0e}")
# Required OSNR of the theory curve
for target in (ber_fec, 1e-3):
    snr_req = brentq(lambda s: ber_16qam(s) - target, 1, 1e4)
    print(
        f"theory: required OSNR {10 * np.log10(snr_req * baud / b_ref):.1f} dB at BER {target:.0e}"
    )
ax.axhline(ber_fec, color="0.5", ls="--", lw=1)
ax.set_xlabel("OSNR in 0.1 nm (dB)")
ax.set_ylabel("BER")
ax.set_ylim(1e-5, 0.5)
ax.legend()
plt.show()

47 GHz front end: required OSNR 23.9 dB at BER 2e-02
47 GHz front end: required OSNR 28.9 dB at BER 1e-03
80 GHz front end: required OSNR 23.6 dB at BER 2e-02
80 GHz front end: required OSNR 28.1 dB at BER 1e-03
theory: required OSNR 22.5 dB at BER 2e-02
theory: required OSNR 26.4 dB at BER 1e-03

# Constellations at the OSNR closest to the FEC threshold, for both front ends
fig, ax = plt.subplots(1, 2, figsize=(9, 4.5))
for a, (front, rows) in zip(ax, results.items()):
    osnr, ber, res = min(
        rows, key=lambda r: abs(np.log10(max(r[1], 1e-9)) - np.log10(ber_fec))
    )
    o = res["X"]["symbols"]
    a.plot(o.real, o.imag, ".", ms=1, alpha=0.2)
    a.set_aspect("equal")
    a.set_xlim(-1.5, 1.5)
    a.set_ylim(-1.5, 1.5)
    a.set_title(f"{labels[front]}\nOSNR {osnr:.1f} dB, BER {ber:.1e}")
plt.show()

At the FEC threshold the two front ends need 1.1 to 1.4 dB more OSNR than the theory curve. The 80 GHz front end gains only 0.3 dB over the 47 GHz front end. At $10^{-3}$ the gap grows to 1.7 and 2.5 dB. There the link has OSNR to spare, and the residual intersymbol interference of the narrower front end matters more. An 800G link of this kind runs at the FEC threshold, so the 47 GHz front end of the receiver page stays.

The remaining penalty comes from the chips and the converters together:

the quantization of the 8-bit DAC and ADC
the TIA noise
the phase noise that the carrier recovery leaves
the imbalance of the hybrids

From an OSNR of 33 dB the 80 GHz front end makes no error in about 99,000 bits.

Summary

Quantity	Value
Line rate	960 Gb/s, DP-16QAM at 120 GBd
Transmitter	LNOI400, 10 mm T-rail modulators, $V_\pi$ 2.95 V, EO bandwidth 101 GHz
Receiver	SOI, SWG edge coupler, Sacher PSR, MMI hybrid, Ge photodiode, 47 GHz front end
Channel	80 km SMF, 1 ps PMD, EDFA with a 5 dB noise figure
Required OSNR at $2 \times 10^{-2}$	23.9 dB (47 GHz), 23.6 dB (80 GHz), theory 22.5 dB
Required OSNR at $10^{-3}$	28.9 dB (47 GHz), 28.1 dB (80 GHz), theory 26.4 dB