Add polarization multiplexing and the blind butterfly equaliser

Two independent 16-QAM tributaries on orthogonal polarizations of one
wavelength: 256 Gb/s at 32 GBd. This is the increment that finally exercises
Ey, which has been in Band since the first commit for exactly this purpose, and
the one that makes the existing PMD model an impairment rather than a curiosity.

The claim worth testing is the failure case. A fibre rotates the launched state
arbitrarily, so a dual-polarization receiver's two branches carry *mixtures*
rather than channels -- and past a few degrees neither is recoverable at all.
Not degraded data: no data. Both halves are asserted, unrecoverable without the
equaliser and error-free with it, at rotations through the 45-degree worst case.

New blocks: Splitter, because a graph edge fans out for free and an optical
field does not, and a dual-pol transmitter really does pay 3 dB;
PolarizationCombiner, which is deliberately the opposite rule to the wavelength
Combiner -- that one refuses two inputs at one frequency because co-located
carriers interfere, this one requires it because orthogonal ones do not;
PolarizationRotator as the impairment; DualPolarizationReceiver;
ButterflyEqualizer.

The equaliser is blind -- no training sequence anywhere in the link. Two stages,
because one does not serve QAM: Godard's single radius opens the eye from
nothing but converges to a compromise no symbol sits on, then radius-directed
error closes it. On QPSK the stages coincide, because QPSK has one radius.

Two things the falsification pass changed rather than confirmed.

A 45-degree rotation mixes the tributaries exactly half and half, which leaves a
symmetric filter initialisation equidistant from both valid solutions -- a
saddle, not a minimum. The adaptation stalls there, and running it longer makes
it *worse*, which is what identified it as a saddle rather than slow
convergence. A small deterministic tilt off the symmetry axis removes it; 64-QAM
at 45 degrees fails without it and lands on the noise floor with it, and a test
guards it.

The other correction is to my own work. The docstring claimed a textbook
singularity guard -- re-seeding the second filter row as the unitary complement
of the first -- was what kept both rows from converging onto the same tributary.
Sabotaging it changed nothing, so it was measured properly across ninety channel
and format combinations: it helped eight times, hurt nine, and did nothing in
the rest. It has been removed and the docstring now says what was measured. A
confident comment over a coin flip is worse than no comment.

Six sabotages, all caught after those two fixes: dropping the tilt, never
leaving the CMA stage, a non-unitary rotator, a combiner that puts both arms on
one axis, a receiver that detects only X, and a filter with no cross terms.

Also honest about what stays unresolved: nothing blind labels the tributaries,
so a channel past 45 degrees separates them perfectly and delivers them the
other way round. A real link recovers the pairing from framing; the tests and
the example recover it by measuring against both references, and say so.

31 components, 450 tests.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

Author: ehsun-sh

README.md

@@ -145,6 +145,31 @@ subtracting a constant or a line, and it puts a ceiling on SNR that no power bud
 search of Pfau et al. Both halves of that claim are [asserted](tests/test_dsp.py) — the second
 would be meaningless without the first.
 
+### Dual polarization
+
+`python examples/dualpol_link.py` puts two independent 16-QAM tributaries on orthogonal
+polarizations of one wavelength — **256 Gb/s** — and rotates the state the way a fibre does:
+
+```
+rotation    without equaliser              with equaliser
+    0 deg  EVM    2.5 /    2.5 %      0 err      EVM 2.50 / 2.57 %    0 err
+   15 deg  EVM   28.1 /   29.2 %   1416 err      EVM 2.51 / 2.58 %    0 err
+   30 deg  EVM  218.0 /  123.5 %   6217 err      EVM 2.52 / 2.51 %    0 err
+   45 deg  EVM  278.5 /  365.2 %   6763 err      EVM 2.48 / 2.55 %    0 err
+   72 deg  EVM   85.7 /  165.6 %   5523 err      EVM 2.54 / 2.56 %    0 err  (swapped)
+   90 deg  EVM    2.5 /    2.5 %      0 err      EVM 2.54 / 2.48 %    0 err  (swapped)
+```
+
+Past a few degrees the unequalised branches are not degraded — they carry no recoverable data at
+all, because each is a *mixture* of both tributaries. The
+[butterfly equaliser](src/oosim/dsp.py) separates them blind, with no training sequence anywhere in
+the link. Read the two end rows together: 90° is a clean swap rather than a mixture, so it needs no
+equaliser at all and simply delivers the tributaries the other way round — which is also why
+nothing blind can label them, and why a real link recovers the pairing from framing.
+
+This is the increment that finally exercises `Ey`, which has been in the signal model since the
+first commit for exactly this purpose.
+
 ```python
 result = sweep(graph, {("laser", "power"): [-24.0, -21.0, -18.0]}, runs=8)
 q = result.metric(analyzer, lambda m: m.q_factor)     # shape (points, runs)
@@ -311,7 +336,7 @@ time window, and results are reproducible.
 | **0 — Foundations** ✅ | Signal model, context, port types, component base, registry, scheduler, `.oosim` project format, sweeps, CI | ~1 month |
 | **1 — MVP: linear link** *(essentially done)* | ✅ PRBS → NRZ → laser → MZM → fiber (α + CD) → PIN → filter → eye/Q/BER, validated end to end. **Python only, no GUI.** | ~2–3 months |
 | **1.5 — Nonlinear & amplified** ✅ | Adaptive-step SSFM, Kerr, EDFA with ASE, OSNR, PMD, APD | ~2 months |
-| **2 — Coherent transceiver** *(single-pol done)* | ✅ Gray-coded M-QAM to 256, IQ modulator with bias and quadrature error, 90° hybrid, balanced detection, blind carrier phase recovery, EVM/MER, constellation diagram, validated against closed-form SER · ⬜ dual polarization, pulse shaping, adaptive equaliser | ~3 months |
+| **2 — Coherent transceiver** ✅ | Gray-coded M-QAM to 256, IQ modulator with bias and quadrature error, 90° hybrid, balanced detection, blind carrier phase recovery, **dual polarization with a blind butterfly equaliser**, EVM/MER, constellation diagram, validated against closed-form SER · ⬜ pulse shaping, differential encoding | ~3 months |
 | **3 — GUI & WDM** | Session server, React Flow graph editor, OSA · DWDM + XPM/FWM crosstalk, 400G/800G references, CuPy back-end | ~6 months |
 | **4 — PIC** | Waveguides, ring resonators, MMI, MZI via integration with an existing S-matrix solver; PDK import | — |
 

examples/dualpol_link.py

@@ -0,0 +1,147 @@
+"""A 256 Gb/s dual-polarization coherent link through a rotated channel.
+
+Two independent 16-QAM tributaries share one wavelength on orthogonal
+polarizations. A fibre rotates the launched state arbitrarily, so what the
+receiver's two branches carry is a *mixture* of both tributaries rather than one
+each — and past a small angle neither is recoverable at all. The butterfly
+equaliser is what separates them again.
+
+Run: ``python examples/dualpol_link.py``
+"""
+
+from __future__ import annotations
+
+from oosim import Graph, SimulationContext
+from oosim.components import (
+    ButterflyEqualizer,
+    CarrierRecovery,
+    ConstellationAnalyzer,
+    CWLaser,
+    DualPolarizationReceiver,
+    IQDriver,
+    IQModulator,
+    IQSampler,
+    PolarizationCombiner,
+    PolarizationRotator,
+    PRBSGenerator,
+    QAMMapper,
+    Splitter,
+)
+
+SYMBOL_RATE = 32e9
+BITS_PER_SYMBOL = 4
+
+
+def build(rotation: float, *, equalize: bool, sequence_length: int = 4096) -> tuple:
+    ctx = SimulationContext(
+        bit_rate=SYMBOL_RATE,
+        samples_per_symbol=4,
+        sequence_length=sequence_length,
+        seed=2026,
+        precision="double",
+    )
+    graph = Graph(ctx)
+    laser = graph.add(CWLaser(power=0.0, linewidth=100.0, label="tx"))
+    splitter = graph.add(Splitter(2, label="sp"))
+    graph.connect(laser, splitter["in"])
+
+    mappers, modulators = {}, {}
+    for index, axis in enumerate(("x", "y")):
+        prbs = graph.add(
+            PRBSGenerator(
+                order=23.0 if axis == "x" else 15.0,
+                bits_per_symbol=float(BITS_PER_SYMBOL),
+                label=f"prbs_{axis}",
+            )
+        )
+        mapper = graph.add(QAMMapper(bits_per_symbol=float(BITS_PER_SYMBOL), label=f"map_{axis}"))
+        driver = graph.add(IQDriver(label=f"drv_{axis}"))
+        modulator = graph.add(IQModulator(label=f"mod_{axis}"))
+        graph.chain(prbs, mapper, driver)
+        graph.connect(splitter[f"out{index}"], modulator["optical_in"])
+        graph.connect(driver["i"], modulator["i"])
+        graph.connect(driver["q"], modulator["q"])
+        mappers[axis], modulators[axis] = mapper, modulator
+
+    combiner = graph.add(PolarizationCombiner(label="pbc"))
+    graph.connect(modulators["x"], combiner["x"])
+    graph.connect(modulators["y"], combiner["y"])
+    rotator = graph.add(PolarizationRotator(angle=rotation, phase=25.0, label="rot"))
+    graph.connect(combiner, rotator["in"])
+
+    lo = graph.add(CWLaser(power=13.0, linewidth=100.0, label="lo"))
+    receiver = graph.add(DualPolarizationReceiver(label="rx"))
+    graph.connect(rotator, receiver["in"])
+    graph.connect(lo, receiver["lo"])
+
+    samplers = {}
+    for axis in ("x", "y"):
+        sampler = graph.add(IQSampler(label=f"smp_{axis}"))
+        graph.connect(receiver[f"{axis}i"], sampler["i"])
+        graph.connect(receiver[f"{axis}q"], sampler["q"])
+        graph.connect(mappers[axis]["out"], sampler["reference"])
+        samplers[axis] = sampler
+
+    if equalize:
+        equalizer = graph.add(ButterflyEqualizer(label="eq"))
+        graph.connect(samplers["x"]["out"], equalizer["x"])
+        graph.connect(samplers["y"]["out"], equalizer["y"])
+        sources = {"x": equalizer["x_out"], "y": equalizer["y_out"]}
+    else:
+        sources = {axis: samplers[axis]["out"] for axis in ("x", "y")}
+
+    analyzers = {}
+    for axis in ("x", "y"):
+        recovery = graph.add(CarrierRecovery(label=f"cr_{axis}"))
+        graph.connect(sources[axis], recovery["in"])
+        for reference in ("x", "y"):
+            analyzer = graph.add(
+                ConstellationAnalyzer(ignore_edges=128.0, label=f"vsa_{axis}{reference}")
+            )
+            graph.connect(recovery["out"], analyzer["in"])
+            graph.connect(mappers[reference]["out"], analyzer["reference"])
+            analyzers[axis + reference] = analyzer
+    return graph, analyzers
+
+
+def measure(rotation: float, *, equalize: bool) -> tuple:
+    graph, analyzers = build(rotation, equalize=equalize)
+    results = graph.run(keep=[])
+    taken = {key: results[a] for key, a in analyzers.items()}
+    # Nothing blind labels the tributaries, so a channel that swaps them is
+    # separated correctly and delivered the other way round. Framing resolves
+    # this in a real link; both references resolve it here.
+    direct = taken["xx"].symbol_errors + taken["yy"].symbol_errors
+    swapped = taken["xy"].symbol_errors + taken["yx"].symbol_errors
+    if direct <= swapped:
+        return taken["xx"], taken["yy"], False
+    return taken["xy"], taken["yx"], True
+
+
+def main() -> None:
+    rate = SYMBOL_RATE * BITS_PER_SYMBOL * 2 / 1e9
+    print(
+        f"Dual-polarization 16-QAM, {SYMBOL_RATE / 1e9:.0f} GBd x 2 pol = {rate:.0f} Gb/s\n"
+        f"Two independent tributaries on one wavelength, through a rotated channel.\n"
+    )
+    print("rotation    without equaliser              with equaliser")
+    print("-" * 70)
+    for rotation in (0.0, 15.0, 30.0, 45.0, 72.0, 90.0):
+        off_x, off_y, _ = measure(rotation, equalize=False)
+        on_x, on_y, swapped = measure(rotation, equalize=True)
+        note = "  (tributaries swapped)" if swapped else ""
+        print(
+            f"{rotation:5.0f} deg  EVM {off_x.evm * 100:6.1f} / {off_y.evm * 100:6.1f} %"
+            f"  {off_x.symbol_errors + off_y.symbol_errors:5} err"
+            f"      EVM {on_x.evm * 100:4.2f} / {on_y.evm * 100:4.2f} %"
+            f"  {on_x.symbol_errors + on_y.symbol_errors:3} err{note}"
+        )
+    print(
+        "\nPast a few degrees the unequalised branches are not degraded — they carry\n"
+        "no recoverable data at all. Separating them is linear algebra the receiver\n"
+        "has to learn blind, with no training sequence anywhere in the link."
+    )
+
+
+if __name__ == "__main__":
+    main()

src/oosim/components/__init__.py

@@ -9,29 +9,43 @@
 
 from .amplifiers import EDFA
 from .analyzers import BERAnalyzer, ConstellationAnalyzer, ConstellationDiagram, EyeDiagram
-from .coherent import CarrierRecovery, CoherentReceiver, IQSampler
+from .coherent import (
+    CarrierRecovery,
+    CoherentReceiver,
+    DualPolarizationReceiver,
+    IQSampler,
+)
 from .detectors import APDPhotodiode, PINPhotodiode
+from .dsp import ButterflyEqualizer
 from .electrical import DCVoltage, IQDriver, NRZDriver, PRBSGenerator
 from .fiber import Fiber
 from .filters import ElectricalFilter
 from .mapping import QAMMapper
 from .meters import OSNRMeter, PowerMeter
 from .modulators import IQModulator, MachZehnderModulator
-from .passive import Attenuator, Combiner
+from .passive import (
+    Attenuator,
+    Combiner,
+    PolarizationCombiner,
+    PolarizationRotator,
+    Splitter,
+)
 from .sources import CWLaser, GaussianPulse, SechPulse
 
 __all__ = [
     "EDFA",
     "APDPhotodiode",
     "Attenuator",
     "BERAnalyzer",
+    "ButterflyEqualizer",
     "CWLaser",
     "CarrierRecovery",
     "CoherentReceiver",
     "Combiner",
     "ConstellationAnalyzer",
     "ConstellationDiagram",
     "DCVoltage",
+    "DualPolarizationReceiver",
     "ElectricalFilter",
     "EyeDiagram",
     "Fiber",
@@ -44,7 +58,10 @@
     "OSNRMeter",
     "PINPhotodiode",
     "PRBSGenerator",
+    "PolarizationCombiner",
+    "PolarizationRotator",
     "PowerMeter",
     "QAMMapper",
     "SechPulse",
+    "Splitter",
 ]

src/oosim/components/coherent.py

@@ -7,6 +7,8 @@
 
 from __future__ import annotations
 
+import math
+
 import numpy as np
 
 from ..component import BoolParam, Component, Param, PortType
@@ -135,6 +137,99 @@ def run(self, ctx: SimulationContext, inputs: dict[str, Signal]) -> dict[str, Si
         }
 
 
+class DualPolarizationReceiver(CoherentReceiver):
+    """Polarization-diverse coherent receiver: four photocurrents, not two.
+
+    The signal is split by a polarization beam splitter, the LO is split with it,
+    and each pair goes to its own 90-degree hybrid. The result is the *complete*
+    optical field: two complex quadratures on each of two orthogonal states, from
+    which any launched polarization can be reconstructed by linear combination.
+
+    That last clause is the point. A single-polarization receiver goes deaf when
+    the fibre rotates the signal onto the orthogonal state; this one cannot,
+    because it measures both states and lets the DSP recombine them. It is why a
+    coherent link can carry two independent tributaries at one wavelength while a
+    direct-detection link carries one.
+
+    The PBS axes are the receiver's own, and they bear no relation to whatever
+    the transmitter launched — a fibre rotates the state arbitrarily and drifts.
+    So ``x`` and ``y`` here are *not* the two tributaries. Separating those is
+    :class:`~oosim.components.dsp.ButterflyEqualizer`'s job, and without it these
+    outputs are two mixtures rather than two channels.
+
+    The noise convention is inherited unchanged: each hybrid sees half the LO
+    power, so the per-quadrature shot variance is half the single-polarization
+    receiver's, and splitting the signal likewise halves what each branch
+    receives. The two effects cancel in the SNR, which is why polarization
+    multiplexing doubles capacity at the same SNR per tributary rather than
+    costing 3 dB.
+    """
+
+    display_name = "Dual-Pol Coherent Receiver"
+    category = "Receivers"
+
+    inputs = {"in": PortType.OPTICAL, "lo": PortType.OPTICAL}
+    outputs = {
+        "xi": PortType.ELECTRICAL,
+        "xq": PortType.ELECTRICAL,
+        "yi": PortType.ELECTRICAL,
+        "yq": PortType.ELECTRICAL,
+    }
+
+    def run(self, ctx: SimulationContext, inputs: dict[str, Signal]) -> dict[str, Signal]:
+        signal: OpticalSignal = inputs["in"]
+        lo: OpticalSignal = inputs["lo"]
+
+        if len(lo.bands) != 1:
+            raise ValueError(
+                f"{self.label}: the local oscillator must be a single band, got "
+                f"{len(lo.bands)}; a coherent receiver mixes against one tone"
+            )
+        reference = lo.bands[0]
+        # The LO is launched at 45 degrees so the PBS delivers equal power to both
+        # hybrids. A real receiver either does this or uses a polarization-
+        # maintaining split; either way each branch gets half.
+        lo_amplitude = math.sqrt(reference.average_power() / 2.0)
+        lo_phase = np.exp(1j * np.angle(reference.Ex.astype(np.complex128)))
+        lo_branch = lo_amplitude * lo_phase
+
+        band = self._nearest_band(signal, reference.f0)
+        responsivity = self.si("responsivity")
+        bandwidth = self.noise_bandwidth(ctx)
+        # Each of the two hybrids sees half the LO.
+        shot_variance = Q_ELECTRON * responsivity * (reference.average_power() / 2.0) * bandwidth
+        thermal_variance = (
+            4.0 * K_BOLTZMANN * self.si("temperature") * bandwidth / self.si("load_resistance")
+            if self.si("load_resistance") > 0.0
+            else 0.0
+        )
+
+        out: dict[str, Signal] = {}
+        for axis, field in (("x", "Ex"), ("y", "Ey")):
+            if band is None:
+                mix = np.zeros(ctx.num_samples, dtype=np.complex128)
+            else:
+                beat = np.exp(2j * np.pi * (band.f0 - reference.f0) * ctx.time_axis())
+                mix = getattr(band, field).astype(np.complex128) * np.conj(lo_branch) * beat
+
+            for quadrature, part in (("i", mix.real), ("q", mix.imag)):
+                current = responsivity * part
+                if self.shot_noise and shot_variance > 0.0:
+                    rng = ctx.rng(type(self).__name__, self.label, "shot", axis, quadrature)
+                    current = current + rng.normal(
+                        0.0, math.sqrt(shot_variance), size=ctx.num_samples
+                    )
+                if self.thermal_noise and thermal_variance > 0.0:
+                    rng = ctx.rng(type(self).__name__, self.label, "thermal", axis, quadrature)
+                    current = current + rng.normal(
+                        0.0, math.sqrt(thermal_variance), size=ctx.num_samples
+                    )
+                out[axis + quadrature] = ElectricalSignal(
+                    samples=current.astype(ctx.real_dtype), fs=ctx.sample_rate, unit="A"
+                )
+        return out
+
+
 class CarrierRecovery(Component):
     """Blind carrier phase recovery, by phase search.