Let the polarizations see each other

The Kerr coupling was scalar per polarization: a channel was modulated by its
neighbours' co-polarized power and by nothing else, so a channel polarized
across it counted for zero. In an isotropic medium orthogonal power counts for
exactly two thirds — which, because co-polarized cross-phase modulation already
carries its factor of two, makes orthogonal cross-phase modulation exactly one
third of co-polarized. Measured on unmodulated carriers over a lossless
dispersionless span, where the phase is a closed form: 2.000 gPL against 0.667.

One coefficient, three consequences. A channel's own orthogonal component
modulates it at two thirds the rate, so power split evenly across the axes turns
each by 0.5 + (2/3)*0.5 instead of 0.5. The two axes stop accumulating the same
phase, so nonlinear birefringence falls by a factor of three. And that
difference is then driven by a *neighbour's* power, which is cross-polarization
modulation — it arrives with the coupling rather than being added beside it.

The solver takes a polarization label per field rather than a second field list,
so the flat argument and the flat return survive and an unlabelled call is the
scalar model term for term. The step-size bound had to learn the same weighting,
or the two paths picked different steps and a single-axis run stopped being
bit-identical to itself; written per field it is now exact rather than an upper
bound, and reduces to the old expression algebraically for co-polarized fields.

The value is the fixed-axis one from the chi(3) tensor, not the Manakov 8/9:
this block applies PMD as a separate element rather than interleaving it, so the
fixed-axis form is what is consistent with the rest of it. The coherent
A_x* A_y**2 term is left out, and a test asserts each axis keeps its power to
twelve digits, so that omission is a number rather than a sentence.

Off by default. With the light on one axis it changes nothing, on the samples.

10 of 10 behaviour-changing sabotages caught. An eleventh — keeping an unused
axis as zeros instead of leaving it absent — changed nothing, correctly: the
docstring claimed that distinction was about exactness and it is about one array
allocation per step, so the docstring was wrong and is now right.

743 tests.

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

Author: ehsun-sh

README.md

@@ -447,6 +447,38 @@ noise bins rendered onto one grid, the way an instrument shows them. Its resolut
 not cosmetic — widening it raises the ASE trace decibel for decibel and leaves a carrier exactly
 where it is, which is the clearest demonstration of why OSNR needs a stated reference bandwidth.
 
+## An orthogonal neighbour is not an absent one
+
+The Kerr coupling used to be scalar per polarization: a channel was modulated by its neighbours'
+co-polarized power and by nothing else, so a channel polarized across it counted for zero. In an
+isotropic medium orthogonal power counts for exactly two thirds — which, since co-polarized
+cross-phase modulation already carries its factor of two, makes orthogonal cross-phase modulation
+exactly **one third** of co-polarized. Set `cross_polarization` on the fibre:
+
+| neighbour | axes uncoupled | axes coupled |
+| :--- | ---: | ---: |
+| co-polarized | 2.000 γPL | 2.000 γPL |
+| orthogonal | 0.000 γPL | 0.667 γPL |
+
+The same coefficient does two more things, because it is the same coefficient. A channel's own
+orthogonal component modulates it at two thirds the rate, so power split evenly between the axes
+turns each of them by `(0.5 + ⅔·0.5)·γPL` instead of `0.5·γPL`. And the two axes, no longer
+accumulating the same phase, rotate the state of polarization as the power moves — nonlinear
+birefringence, which falls by a factor of three, and cross-polarization modulation, which is that
+rotation being driven by a *different* channel's power.
+
+The value is the fixed-axis one, from the χ⁽³⁾ tensor rather than from any averaging. A fibre whose
+birefringence scrambles faster than the nonlinearity acts is the Manakov regime instead, where the
+distinction washes into a single 8/9 on the total power; this block applies PMD as a separate
+element rather than interleaving it, so the fixed-axis form is the one consistent with the rest of
+it. The coherent `A_x* A_y²` term, which would exchange power between the axes rather than only
+dephase them, is left out — it is the part that averages away first — and the tests assert the
+axes' powers are unchanged to twelve digits, so that omission is a number rather than a sentence.
+
+Off by default, and with all the light on one axis it changes nothing at all — on the samples,
+which is what makes it safe to leave on for a dual-polarization link and pointless for a
+single-polarization one.
+
 ## D is not one number
 
 Standard fibre gains 0.058 ps/nm/km of dispersion for every nanometre up the

examples/ui_data.json

@@ -724,6 +724,11 @@
      "default": true,
      "doc": "Couple the bands: each is phase-modulated by the others' power"
     },
+    "cross_polarization": {
+     "type": "bool",
+     "default": false,
+     "doc": "Couple the two polarizations: orthogonal power modulates at two thirds"
+    },
     "max_walkoff_slip": {
      "type": "float",
      "default": 0.5,

examples/wdm_nonlinear.py

@@ -67,7 +67,8 @@
     fwm_phase_mismatch,
     walkoff_from_dispersion,
 )
-from maiman.units import C_LIGHT, w_to_dbm, wavelength_to_frequency
+from maiman.signals import Band, OpticalSignal
+from maiman.units import C_LIGHT, dbm_to_w, w_to_dbm, wavelength_to_frequency
 
 ANCHOR = 1550.0  # nm — channel 0
 SPACING = 100e9  # Hz
@@ -308,6 +309,51 @@ def main() -> None:
     print("     where no filter downstream can reach it. That is what a uniform grid")
     print("     costs, and why unequal channel spacing was once a serious proposal.")
 
+    # ------------------------------------------------------------------
+    print("\n  6. Polarization: an orthogonal neighbour is not an absent one")
+    print("     Two unmodulated channels, 100 km, no loss, no dispersion.")
+    print("     The phase channel 0 picks up beyond its own self-phase modulation:\n")
+    ctx = SimulationContext(bit_rate=10e9, samples_per_symbol=16, sequence_length=64, seed=3)
+    power = dbm_to_w(0.0)
+    unit = GAMMA * 1e-3 * power * 100e3  # gamma * P * L, radians
+
+    def carrier(px: float, py: float, index: int) -> Band:
+        samples = ctx.num_samples
+        return Band(
+            Ex=np.full(samples, np.sqrt(px), dtype=np.complex128),
+            Ey=np.full(samples, np.sqrt(py), dtype=np.complex128),
+            f0=wavelength_to_frequency(channel_wavelength(index) * 1e-9),
+            fs=ctx.sample_rate,
+        )
+
+    def turned(bands: tuple[Band, ...], *, coupled: bool) -> float:
+        span = Fiber(
+            length=100.0,
+            attenuation=0.0,
+            dispersion=0.0,
+            nonlinearity=GAMMA,
+            four_wave_mixing=False,
+            cross_polarization=coupled,
+            label="fib",
+        )
+        out = span.run(ctx, {"in": OpticalSignal(bands=bands, noise=())})["out"]
+        return float(np.angle(out.bands[0].Ex[0])) / unit
+
+    print(f"     {'neighbour':>22}  {'axes uncoupled':>15}  {'axes coupled':>13}")
+    print("     " + "-" * 54)
+    for label, neighbour in (
+        ("co-polarized", carrier(power, 0.0, 1)),
+        ("orthogonal", carrier(0.0, power, 1)),
+    ):
+        alone_off = turned((carrier(power, 0.0, 0),), coupled=False)
+        alone_on = turned((carrier(power, 0.0, 0),), coupled=True)
+        off = turned((carrier(power, 0.0, 0), neighbour), coupled=False) - alone_off
+        on = turned((carrier(power, 0.0, 0), neighbour), coupled=True) - alone_on
+        print(f"     {label:>22}  {off:12.3f} gPL  {on:10.3f} gPL")
+    print("     Two and two thirds: orthogonal cross-phase modulation is exactly a")
+    print("     third of co-polarized. Uncoupled it is zero, which is the model")
+    print("     saying a channel across the polarization is not there at all.")
+
 
 if __name__ == "__main__":
     main()

src/maiman/components/fiber.py

@@ -102,14 +102,23 @@ class Fiber(Component):
     the run's generator, which makes a single run reproducible and a sweep with
     repeats an exploration of the distribution.
 
-    Not yet modelled: Raman scattering. The
-    Kerr coupling is scalar per polarization throughout, so a channel is
-    modulated by its neighbours' co-polarized power and not by their orthogonal
-    power, which really contributes a third as much; cross-polarization
-    modulation is absent for the same reason. Mixing products accumulate from
-    span to span in power rather than in field, which understates the coherent
-    build-up a dispersion-managed link produces. Pump depletion is not modelled,
-    which matters only at powers no link is operated at.
+    **The polarizations couple only if asked.** By default each is propagated as
+    its own scalar problem, which is what every result in this project was taken
+    with. Set ``cross_polarization`` and orthogonal power enters the Kerr term at
+    two thirds the co-polarized weight — so a neighbour polarized across the
+    channel modulates it at exactly one third of the rate a co-polarized one
+    does, and the two axes of one channel, no longer accumulating the same phase,
+    rotate the state of polarization as the power moves. With all the light on
+    one axis the setting changes nothing, which is why it is safe to leave on for
+    a dual-polarization link and pointless for a single-polarization one.
+
+    Not yet modelled: Raman scattering. The coherent polarization term
+    ``A_x* A_y**2``, which exchanges power between the axes rather than only
+    dephasing them, is left out — it is the part that averages away first under
+    real birefringence. Mixing products accumulate from span to span in power
+    rather than in field, which understates the coherent build-up a
+    dispersion-managed link produces. Pump depletion is not modelled, which
+    matters only at powers no link is operated at.
     """
 
     display_name = "Optical Fiber"
@@ -140,6 +149,9 @@ class Fiber(Component):
     cross_phase_modulation = BoolParam(
         True, doc="Couple the bands: each is phase-modulated by the others' power"
     )
+    cross_polarization = BoolParam(
+        False, doc="Couple the two polarizations: orthogonal power modulates at two thirds"
+    )
     max_walkoff_slip = Param(
         0.5,
         unit="",
@@ -316,9 +328,16 @@ def _propagate_kerr(
     ) -> tuple[list[tuple[np.ndarray, np.ndarray]], PropagationDiagnostics]:
         """Split-step propagation, with the bands coupled unless told otherwise.
 
-        The two polarizations are propagated as two separate coupled systems,
-        which is what makes the model scalar per polarization: a band's X field
-        is modulated by every other band's X power and by nothing else.
+        By default the two polarizations are propagated as two separate coupled
+        systems, which is what makes the model scalar per polarization: a band's
+        X field is modulated by every other band's X power and by nothing else.
+        Set ``cross_polarization`` and both axes go into one call, labelled, so
+        that orthogonal power enters at :data:`maiman.kernels.ORTHOGONAL_KERR_WEIGHT`.
+
+        The two axes of one band share a linear operator — the same beta2, the
+        same slope, the same walk-off — because birefringence is this block's
+        other business and is applied as a separate element afterwards. What the
+        coupling adds is entirely in the nonlinear step.
         """
         bands = signal.bands
         coupled = self.cross_phase_modulation and len(bands) > 1
@@ -340,26 +359,74 @@ def _propagate_kerr(
             beta2 = [self.beta2_at(band.wavelength) for band in group]
             beta3 = [self.beta3_at(band.wavelength) for band in group]
             walkoff = [self.walkoff_of(band, reference) for band in group]
+
+            if self.cross_polarization:
+                # One system, both axes, labelled. The per-field lists are
+                # doubled rather than special-cased in the solver: X and Y of the
+                # same band see the same linear operator.
+                out, diagnostics = self._solve(
+                    [band.Ex for band in group] + [band.Ey for band in group],
+                    reference.fs,
+                    beta2=beta2 * 2,
+                    beta3=beta3 * 2,
+                    walkoff=walkoff * 2,
+                    polarization=[0] * len(group) + [1] * len(group),
+                    gamma=gamma,
+                    alpha=alpha,
+                    distance=distance,
+                    best=diagnostics,
+                )
+                fields.extend(zip(out[: len(group)], out[len(group) :], strict=True))
+                continue
+
             solved = []
             for axis in ("Ex", "Ey"):
-                out, diag = propagate_coupled_ssfm(
+                out, diagnostics = self._solve(
                     [getattr(band, axis) for band in group],
                     reference.fs,
                     beta2=beta2,
+                    beta3=beta3,
                     walkoff=walkoff,
+                    polarization=None,
                     gamma=gamma,
-                    beta3=beta3,
                     alpha=alpha,
                     distance=distance,
-                    max_nonlinear_phase=self.max_nonlinear_phase,
-                    max_walkoff_slip=self.max_walkoff_slip,
+                    best=diagnostics,
                 )
                 solved.append(out)
-                if diag.steps > diagnostics.steps:
-                    diagnostics = diag
             fields.extend(zip(solved[0], solved[1], strict=True))
         return fields, diagnostics
 
+    def _solve(
+        self,
+        fields: list[np.ndarray],
+        sample_rate: float,
+        *,
+        beta2: list[float],
+        beta3: list[float],
+        walkoff: list[float],
+        polarization: list[int] | None,
+        gamma: float,
+        alpha: float,
+        distance: float,
+        best: PropagationDiagnostics,
+    ) -> tuple[list[np.ndarray], PropagationDiagnostics]:
+        """One call into the solver, keeping whichever diagnostics took more steps."""
+        out, diag = propagate_coupled_ssfm(
+            fields,
+            sample_rate,
+            beta2=beta2,
+            walkoff=walkoff,
+            gamma=gamma,
+            beta3=beta3,
+            polarization=polarization,
+            alpha=alpha,
+            distance=distance,
+            max_nonlinear_phase=self.max_nonlinear_phase,
+            max_walkoff_slip=self.max_walkoff_slip,
+        )
+        return out, diag if diag.steps > best.steps else best
+
     # -- four-wave mixing -------------------------------------------------
 
     def _mix(