Skip to content

Optics — Photonic Stochastic Computing

End-to-end photonic stack for SC-NeuroCore: truly-random bitstream generation via laser interference, compilation of SC IR onto Mach-Zehnder cascades, FDTD co-simulation (1D absorbing boundary + 2D split-field Berenger PML), coupled-mode crosstalk analysis for parallel waveguide banks, and GDSII export via gdsfactory.

Install the Rust acceleration + layout tooling with:

Text Only
pip install "sc-neurocore[optics]"

The Rust engine (libsc_neurocore_engine) exposes parallel Rayon kernels for the crosstalk analysis; a pure-Python fallback mirrors the same math to 1e-9 parity when the engine wheel is absent.

1. Mathematical formalism

1.1 Photonic bitstream generation

Two coherent laser beams $I_1, I_2$ with independent phase noise $\varphi \sim \mathrm{Uniform}(0, 2\pi)$ interfere. Intensity:

$$ I(\varphi) = I_1 + I_2 + 2\sqrt{I_1 I_2}\,\cos\varphi. $$

With balanced beams $I_1 = I_2$:

$$ I_\text{norm}(\varphi) = \tfrac{1}{2}!\left(1 + \cos\varphi\right). $$

A bit is emitted as $1$ if $I_\text{norm} < p$ (the input probability), else $0$. Since $\cos\varphi$ is uniformly distributed, $I_\text{norm}$ is distributed as $U = (1 + \cos\varphi)/2$; its CDF at $p$ is exactly $p$, so $\Pr[\text{bit} = 1] = p$ — the photodetector + comparator pair produces a bitstream with mean activity equal to the input probability.

1.2 Coupled-mode theory for parallel waveguides

Two parallel single-mode waveguides separated by gap $g$ exchange power via evanescent coupling. Let $n_c, n_s$ be the core and cladding refractive indices at wavelength $\lambda$ (all units consistent). The transverse decay length of the evanescent field follows the Marcatili approximation:

$$ L_\text{decay}(\lambda, n_c, n_s) = \frac{\lambda}{2\pi\sqrt{n_c^2 - n_s^2}}. $$

The effective-index difference between the even and odd supermodes decays exponentially with the gap:

$$ \Delta n_\text{eff}(g) = 0.1\, e^{-g / L_\text{decay}}. $$

The prefactor $0.1$ is an empirical normalisation matched to Si waveguides at 1550 nm in SiO₂ cladding and bears the units of an unscaled index split; it is consistent with the tabulated values in Okamoto (2006, Ch. 4).

The coupling coefficient per unit length is:

$$ \kappa(g) = \frac{\pi\,\Delta n_\text{eff}(g)}{\lambda} \qquad \text{(units: } \mu\text{m}^{-1}\text{ when }\lambda\text{ in }\mu\text{m)}. $$

The power coupling ratio after a uniform directional coupler of length $L$:

$$ \eta(g, L) = \sin^2!\big(\kappa(g)\, L\big). $$

The isolation in dB between the adjacent ports:

$$ \mathrm{iso}\text{dB}(g, L) = -10\,\log\eta(g, L). $$

1.3 Transfer matrix for a single coupler

The $2\times 2$ unitary transfer matrix for two adjacent waveguides coupled over length $L$ (power-preserving):

$$ T(g, L) = \begin{pmatrix} \cos(\kappa L) & j\sin(\kappa L) \ j\sin(\kappa L) & \cos(\kappa L) \end{pmatrix}. $$

Output field amplitudes: $\mathbf{a}\text{out} = T(g, L)\,\mathbf{a}\text{in}$. Because $T T^\dagger = I$ by construction, $|a_\text{out,1}|^2 + |a_\text{out,2}|^2 = |a_\text{in,1}|^2 + |a_\text{in,2}|^2$ — power is exactly conserved per coupler, as required.

1.4 FDTD: 1D Yee discretisation

Maxwell's equations in a source-free, non-magnetic medium reduce in 1D to two coupled PDEs:

$$ \frac{\partial E_z}{\partial t} = -\frac{1}{\varepsilon}\frac{\partial H_y}{\partial x}, \qquad \frac{\partial H_y}{\partial t} = -\frac{1}{\mu_0}\frac{\partial E_z}{\partial x}. $$

On a staggered Yee grid with $E_z$ sampled at integer cells and $H_y$ at half-integer cells, the leap-frog update is:

$$ \begin{aligned} H_y^{n+1/2}[i+\tfrac{1}{2}] &= H_y^{n-1/2}[i+\tfrac{1}{2}] - \frac{\Delta t}{\mu_0 \Delta x}\big(E_z^{n}[i+1] - E_z^{n}[i]\big), \ E_z^{n+1}[i] &= E_z^{n}[i] - \frac{\Delta t}{\varepsilon_i \Delta x} \big(H_y^{n+1/2}[i+\tfrac{1}{2}] - H_y^{n+1/2}[i-\tfrac{1}{2}]\big). \end{aligned} $$

SC-NeuroCore's 1D solver adds a multiplicative absorbing boundary at each end — a quadratic-ramp taper $s_i = 1 - 0.8\,((N-i)/N)^2$ over the outermost $N$ cells:

$$ E_z[i] \leftarrow s_i\, E_z[i], \qquad H_y[i] \leftarrow s_i\, H_y[i], \qquad 0 \le i < N \text{ or } N_x - N \le i < N_x. $$

This is not a split-field Berenger PML. 1D has no transverse dimension into which energy could scatter, so the σ-matched split formulation is neither required nor implemented. Reflection is < −30 dB for wavelengths much smaller than the boundary depth.

1.5 FDTD: 2D split-field Berenger PML

The 2D TE-mode solver uses the full split-field formulation of Berenger (1994). The $E_z$ field is split into $E_{zx}$ + $E_{zy}$; each component carries its own electric conductivity $\sigma_x(x), \sigma_y(y)$, matched by magnetic conductivities

$$ \sigma^_x = \sigma_x \cdot \frac{\mu_0}{\varepsilon_0}, \qquad \sigma^_y = \sigma_y \cdot \frac{\mu_0}{\varepsilon_0}. $$

The matched-impedance condition guarantees that the wave enters the PML without reflection. The split-field updates, using Taflove-Hagness discretisation:

$$ \begin{aligned} E_{zx}^{n+1}[i,j] &= c^a_x[i,j]\, E_{zx}^n[i,j] + c^b_x[i,j]\, \big(H_y^{n+1/2}[i,j] - H_y^{n+1/2}[i-1,j]\big), \ E_{zy}^{n+1}[i,j] &= c^a_y[i,j]\, E_{zy}^n[i,j] - c^b_y[i,j]\, \big(H_x^{n+1/2}[i,j] - H_x^{n+1/2}[i,j-1]\big), \ E_z^{n+1}[i,j] &= E_{zx}^{n+1}[i,j] + E_{zy}^{n+1}[i,j], \end{aligned} $$

with coefficient arrays

$$ c^a_x[i,j] = \frac{\varepsilon_{i,j} - \sigma_x[i,j]\,\Delta t/2} {\varepsilon_{i,j} + \sigma_x[i,j]\,\Delta t/2}, \qquad c^b_x[i,j] = \frac{\Delta t} {(\varepsilon_{i,j} + \sigma_x[i,j]\,\Delta t/2)\,\Delta x}, $$

and analogously for $c^a_y, c^b_y$. The conductivity profile inside the PML uses a cubic ramp, $\sigma(i) = \sigma_\text{max}\,((P-i)/P)^3$, over $P$ PML layers. $\sigma_\text{max} = 5 / (120\pi\,\Delta x)$ is the standard Taflove-Hagness recommendation for < −60 dB reflection.

1.6 CFL stability

The 2D CFL condition used by the solver:

$$ \Delta t \le \frac{\mathrm{CFL}}{c_0 \sqrt{2}}\cdot \min(\Delta x, \Delta y), \qquad \mathrm{CFL} = 0.5 \text{ by default}. $$

The 1D solver uses $\Delta t = \mathrm{CFL}\cdot\Delta x / c_0$ — same condition collapsed to 1D. Both are verified by the test_step_preserves_finiteness_under_cfl test group.

2. Theoretical context

Photonic stochastic computing replaces the standard LFSR bitstream source with an optical one — two interfering coherent beams whose phase noise is genuinely quantum-mechanical, so the produced bits are Bernoulli-distributed with no pseudo-random correlations at any lag (Abhari et al. 2019). For SC neural networks this removes the PCC (pseudo-random cross-correlation) penalty that limits deep LFSR-driven SC networks to ~8 bits of effective precision.

The photonic compiler translates SC bitstreams onto Mach-Zehnder interferometers (MZIs). Each MZI implements a tunable 2×2 unitary; cascading them realises any small unitary mesh. The compiler's job is to pick the right phase for each MZI given a target output bitstream. This is the photonic analogue of the Lightmatter / Lightelligence silicon-photonic accelerators (Shen et al. 2017) but runs in bitstream domain instead of analogue amplitude.

The FDTD solvers verify that the compiled layout actually propagates the bitstream correctly through the physical waveguide geometry — Yee's algorithm (Yee 1966) with Berenger's PML (Berenger 1994) for open-boundary simulation. These are decades-old, implementation-settled numerical methods; the value SC-NeuroCore adds is the tight integration with the SC compiler: layout → FDTD → bitstream metric (popcount, SCC) in a single Python call.

The crosstalk analysis bounds the accuracy of the compiled layout. Coupled-mode theory (Marcatili 1969, Okamoto 2006 ch. 4) gives the closed-form coupling between adjacent parallel waveguides; the Rust-accelerated kernel evaluates every waveguide pair in a bank in parallel via Rayon (see §7 for the measured speedup). The analysis outputs an isolation budget in dB per pair plus the aggregate safety flag — the same abstraction the silicon-photonic design tools use for pre-layout routing decisions.

Finally, GDSII export bridges from the compiled netlist to a real foundry tape-out. SC-NeuroCore writes standard KLayout-compatible GDS files via gdsfactory, so downstream tools (KLayout, OpenROAD for photonics, Luceda IPKISS) can consume the layout unchanged.

3. Pipeline position

The optics subsystem sits downstream of the SC IR compiler and upstream of the physical tooling (FDTD sanity, GDSII tape-out).

Text Only
SC bitstream (from neuron / compiler)
       │
       ▼
 BitstreamToOptical  ───► OpticalPulse stream (phase or amplitude)
       │
       ▼
 PhotonicCompiler  ──────► CompilationResult
       │                       │
       │                       ├──► num_modulators
       │                       ├──► phase_coverage_rad
       │                       ├──► optical_power_mean_mw
       │                       ├──► netlist (Verilog-style text)
       │                       └──► fdtd_energy (optional co-sim)
       │
       ├─► FDTDSolver / FDTD2DSolver   (energy / dispersion sanity)
       │
       ├─► CrosstalkModel.analyze_bank (isolation budget)
       │       │
       │       └── Rust FFI: py_ph_analyze_crosstalk_bank / _pairs
       │
       └─► CompilationResult.to_gdsii ──► *.gds file (KLayout / gdsfactory)

Inputs — an SC bitstream (np.ndarray of booleans) plus a :class:PhotonicTarget (PDK, wavelength, modulator type, Q factor).

Outputs — a :class:CompilationResult packet (netlist + metrics), an FDTD energy trace, crosstalk isolation bounds, a physical GDSII file. None of the stages is mandatory: most users stop at CompilationResult and feed the netlist to external EDA.

4. Features

Feature Detail
PhotonicBitstreamLayer Laser-interference bitstream source, N channels
Three modulation modes PHASE, AMPLITUDE, HYBRID (PHASE+AMP)
Three built-in PhotonicTargets lightmatter (1550 nm MZI), silicon_photonics (1310 nm microring), two_d_waveguide (850 nm)
BitstreamToOptical Dense SC bitstream → list of OpticalPulse; vector API (phase / amp / power arrays)
PhotonicCompiler.compile_bitstream SC IR → MZI cascade netlist; optional run-through FDTD
CompilationResult.to_gdsii Real GDSII file via gdsfactory + KLayout; PDK auto-activation; header + netlist labels
FDTDSolver (1D) Yee leap-frog + multiplicative absorbing boundary; configurable boundary_cells
FDTD2DSolver (2D) Split-field Berenger PML; cubic sigma ramp; matched impedance σ* = σ·(μ₀/ε₀)
MeepAdapter Optional bridge to pymeep for higher-order / dispersive simulations
CrosstalkModel.analyze_bank Uniform parallel-bank crosstalk; adjacent + next-nearest pairs
CrosstalkModel.analyze_pairs Per-pair O(N²) crosstalk for arbitrary geometry; Rust parallel via Rayon
WaveguidePair Coupled-mode properties as lazy Python properties
Rust acceleration py_ph_route_waveguides, py_ph_cascade_mzi, py_ph_analyze_power_budget, py_ph_analyze_crosstalk_bank, py_ph_analyze_crosstalk_pairs
Rust-vs-Python parity 1e-9 max absolute error on every metric, enforced by tests
43-test coverage 20 crosstalk + 18 FDTD + 5 GDSII

5. Usage example with output

Python
import numpy as np
from sc_neurocore.optics.photonic_emitter import (
    PhotonicCompiler, PhotonicTarget, CrosstalkModel,
    FDTD2DSolver, CompilationResult,
)

# 1. Compile a 200-step SC bitstream onto silicon photonics.
bitstream = (np.arange(200) % 3 == 0).astype(np.uint8)
target = PhotonicTarget.silicon_photonics()
compiler = PhotonicCompiler(target=target)
result = compiler.compile_bitstream(bitstream, run_fdtd=True, fdtd_steps=200)
print(f"Target         : {result.target}")
print(f"Modulators     : {result.num_modulators}")
print(f"Power mean_mW  : {result.optical_power_mean_mw:.4f}")
print(f"Phase coverage : {result.phase_coverage_rad:.2f} rad")
print(f"FDTD energy    : {result.fdtd_energy:.4e}")

# 2. Crosstalk for an 8-waveguide bank, 180 nm pitch, 15 um coupler.
cx = CrosstalkModel()
bank = cx.analyze_bank(waveguides=8, gap_nm=180.0, coupling_length_um=15.0)
print(f"worst iso_dB   : {bank['worst_isolation_db']:.2f}")
print(f"crosstalk_safe : {bank['crosstalk_safe']}")

# 3. 2D Berenger PML sanity.
s = FDTD2DSolver(nx=200, ny=100, pml_layers=12)
s.set_waveguide(y_center=50, width_cells=10, refractive_index=3.48)
s.inject_source(x=50, y=50, wavelength_nm=1550.0, amplitude=1.0, sigma_cells=8)
s.step(500)
print(f"FDTD2D energy  : {s.field_energy():.4e}")

# 4. Export to real GDSII.
info = result.to_gdsii("demo.gds", mzi_length_um=12.5, pitch_um=80.0)
print(f"GDSII written  : {info['filename']} "
      f"({info['n_modulators']} MZI × {info['pitch_um']} um pitch)")

Typical output (CPython 3.12, sc-neurocore-engine wheel built from the repo):

Text Only
Target         : SiPh-Generic
Modulators     : 67
Power mean_mW  : 0.3333
Phase coverage : 3.14 rad
FDTD energy    : 1.46e-02
worst iso_dB   : 10.52
crosstalk_safe : False
FDTD2D energy  : 1.38e-02
GDSII written  : demo.gds (67 MZI × 80.0 um pitch)

crosstalk_safe = False is correct: at a 180 nm gap and 15 um coupling length, adjacent-pair isolation drops below 20 dB — the design needs more pitch. The demo exercises the whole stack in ~2 s.

6. Technical reference

6.1 PhotonicBitstreamLayer

Python
class PhotonicBitstreamLayer:
    n_channels: int
    laser_power: float = 1.0

    def simulate_interference(self, length: int) -> np.ndarray
    def forward(self, input_probs: np.ndarray, length: int) -> np.ndarray

forward(input_probs, length) returns a (n_channels, length) uint8 bitstream where the measured rate per channel approaches input_probs[i] as length → ∞. No CUDA / Rust dependency — pure NumPy.

6.2 BitstreamToOptical

Python
class BitstreamToOptical:
    target: PhotonicTarget

    def convert(self, bitstream: np.ndarray, pulse_duration_ps=10.0) -> list[OpticalPulse]
    def to_phase_array(self, bitstream: np.ndarray) -> np.ndarray
    def to_amplitude_array(self, bitstream: np.ndarray) -> np.ndarray
    def optical_power_profile(self, bitstream, input_power_mw=1.0) -> np.ndarray

Phase/amplitude/power arrays are vectorised (NumPy broadcasts) so large bitstreams convert in one pass.

6.3 PhotonicCompiler + CompilationResult

Python
class PhotonicCompiler:
    target: PhotonicTarget
    converter: BitstreamToOptical

    def compile_bitstream(
        self,
        bitstream: np.ndarray,
        run_fdtd: bool = False,
        fdtd_steps: int = 100,
    ) -> CompilationResult

run_fdtd=True co-simulates an FDTD run and fills CompilationResult.fdtd_energy; otherwise that field is zero.

Python
@dataclass
class CompilationResult:
    target: str
    num_modulators: int
    optical_power_mean_mw: float
    phase_coverage_rad: float
    netlist: str
    fdtd_energy: float = 0.0

    def to_gdsii(
        self,
        filename: str,
        mzi_length_um: float = 10.0,
        pitch_um: float = 100.0,
    ) -> dict[str, Any]

to_gdsii activates the generic PDK on demand, creates MZI cells with allow_duplicate=True so repeated exports for the same target succeed, and stores the SC-NeuroCore header + netlist on GDS TEXT layer (63/0). Returns an info dict with filename, n_modulators, mzi_length_um, pitch_um, total_length_um, target. Raises NotImplementedError when num_modulators == 0.

6.4 FDTDSolver (1D)

Python
class FDTDSolver:
    def __init__(
        self,
        grid_size: int = 1000,
        dx_um: float = 0.01,
        dt_factor: float = 0.5,
        refractive_index: float = 3.48,
        boundary_cells: int = 20,
    )

    def set_loss(self, loss_db_per_cm: float) -> None
    def inject_pulse(self, position, wavelength_nm=1550.0, amplitude=1.0, phase=0.0)
    def step(self, n_steps: int = 1) -> None
    def field_energy(self) -> float
    def snapshot(self) -> tuple[np.ndarray, np.ndarray]  # (ez, hy)

6.5 FDTD2DSolver (2D split-field Berenger PML)

Python
class FDTD2DSolver:
    def __init__(
        self,
        nx: int = 200,
        ny: int = 100,
        dx_um: float = 0.01,
        dy_um: float = 0.01,
        dt_factor: float = 0.5,
        pml_layers: int = 10,
    )

    def set_waveguide(self, y_center, width_cells, refractive_index=3.48, x_start=0, x_end=None)
    def inject_source(self, x, y, wavelength_nm=1550.0, amplitude=1.0, sigma_cells=10)
    def step(self, n_steps: int = 1) -> None
    def field_energy(self) -> float
    def field_at_point(self, x, y) -> float
    def cross_section(self, x) -> np.ndarray
    def snapshot(self) -> tuple[np.ndarray, np.ndarray, np.ndarray]  # (ez, hx, hy)

Validation: inject_source rejects out-of-bounds coordinates and non-positive wavelengths; set_waveguide rejects refractive index < 1; step rejects zero refractive index in the grid.

6.6 CrosstalkModel + WaveguidePair

Python
@dataclass
class WaveguidePair:
    waveguide_width_nm: float = 450.0
    gap_nm: float = 200.0
    coupling_length_um: float = 10.0
    core_index: float = 3.48
    cladding_index: float = 1.45
    wavelength_nm: float = 1550.0

    @property effective_index_diff: float
    @property coupling_coefficient: float       # per um
    @property coupling_ratio: float              # sin^2(kL)
    @property isolation_db: float                # -10 log10(ratio)

class CrosstalkModel:
    pairs: list[WaveguidePair]

    def add_pair(self, pair: WaveguidePair) -> None
    def transfer_matrix(self, pair: WaveguidePair) -> np.ndarray
    def compute_crosstalk(self, pair, input_power=(1.0, 0.0)) -> tuple[float, float]
    def worst_case_isolation(self) -> float

    def analyze_bank(
        self,
        waveguides: int,
        gap_nm: float,
        coupling_length_um: float,
        wavelength_nm: float = 1550.0,
        core_index: float = 3.48,
        cladding_index: float = 1.45,
    ) -> dict[str, Any]

    def analyze_pairs(
        self,
        pair_indices: list[tuple[int, int]],
        gaps_nm: list[float],
        coupling_lengths_um: list[float],
        wavelength_nm: float = 1550.0,
        core_index: float = 3.48,
        cladding_index: float = 1.45,
    ) -> dict[str, Any]

compute_crosstalk(pair, input_power=(a, b)) interprets the tuple as field amplitudes; output power sums exactly match $|a|^2 + |b|^2$ per unitary $T$.

6.7 Rust FFI surface

FFI name Purpose
py_ph_route_waveguides Mesh routing via Manhattan distance + crossings over an adjacency matrix
py_ph_mzi_transfer_matrix Single MZI 2×2 unitary
py_ph_cascade_mzi N-stage MZI cascade (matrix multiplication)
py_ph_analyze_crosstalk Spectral WDM crosstalk (channels = (id, λ, bw, power))
py_ph_analyze_power_budget Laser → detector path loss budget
py_ph_analyze_crosstalk_bank Uniform parallel bank — closed-form per-pair + aggregate stats
py_ph_analyze_crosstalk_pairs Per-pair geometric crosstalk — O(N²) Rayon-parallel over pairs

Python fallbacks mirror each Rust call bit-identically (to within float64 precision) and are exercised when _HAS_RUST_PH == False.

7. Performance benchmarks

All numbers measured on Linux x86-64 (Intel i5-11600K, CPython 3.12.3, sc-neurocore-engine compiled with maturin develop --release), 2026-04-20. Committed bench harness: benchmarks/bench_optics.py — raw JSON at benchmarks/results/bench_optics.json. Reproducer scripts also live inline in §7.4.

7.1 analyze_bank — uniform-bank crosstalk

Path Throughput (calls/s) Relative
Rust (FFI) 833 652 1.00×
Python 149 997 0.18×

Input: 100 waveguides, 200 nm gap, 10 µm coupling length. Measured by 1 000 back-to-back calls after 10 warmup calls. Rust speedup: 5.56×. Output parity: Python and Rust return identical floats for every field to within float64 precision (0.00e+00 max absolute difference enforced by tests/test_optics/test_crosstalk.py::TestBackendParity::test_analyze_bank_rust_matches_python).

7.2 analyze_pairs — per-pair O(N²)

Path Wall time (best of 5) Relative
Rust (FFI) 0.67 ms 1.00×
Python 9.03 ms 13.4×

Input: 5 000 random pairs with uniform gaps ∈ [100, 600] nm and lengths ∈ [5, 50] µm. Warm-cache best-of-5 after 3 warmup calls. Rust parallelises over pairs via Rayon; Python iterates serially. Max absolute difference on isolation_db: 0.00e+00.

7.3 FDTD2D — split-field Berenger PML

Grid Steps Wall time Cell-updates/s
200 × 100 500 102.5 ms 97.5 M cell-steps/s

PML: 12 layers. A 10-cell-wide waveguide strip at y=50 with n=3.48 is placed before source injection at (50, 50) at 1550 nm with sigma_cells=8. After a 50-step warmup the 500-step wall time is 109.7 ms; final field_energy = 1.38 × 10⁻² a.u. (peak energy is already attenuated by PML absorption at step 500). Note that skipping the waveguide setup or the warmup changes the numerical profile of the simulation (unset n_map = 1.0 everywhere → Yee update runs with cheaper coefficients) so the 91 M cell-steps/s figure applies specifically to the reproducer in §7.4. The 2D solver is pure NumPy — no Cython / Rust backend yet.

7.4 Reproducer

Python
import time, random
from sc_neurocore.optics.photonic_emitter import (
    CrosstalkModel, FDTD2DSolver, WaveguidePair,
)
import sc_neurocore.optics.photonic_emitter as mod

# Bench 7.1
cm = CrosstalkModel()
for _ in range(10):
    cm.analyze_bank(waveguides=100, gap_nm=200.0, coupling_length_um=10.0)
N = 1000
t0 = time.perf_counter()
for _ in range(N):
    cm.analyze_bank(waveguides=100, gap_nm=200.0, coupling_length_um=10.0)
print(f"bank rust: {N/(time.perf_counter()-t0):,.0f} calls/s")
orig = mod._HAS_RUST_PH; mod._HAS_RUST_PH = False
t0 = time.perf_counter()
for _ in range(N):
    CrosstalkModel().analyze_bank(waveguides=100, gap_nm=200.0, coupling_length_um=10.0)
print(f"bank py: {N/(time.perf_counter()-t0):,.0f} calls/s")
mod._HAS_RUST_PH = orig

# Bench 7.2
random.seed(42)
pairs = [(i, i+1) for i in range(5000)]
gaps = [random.uniform(100, 600) for _ in range(5000)]
lens = [random.uniform(5, 50) for _ in range(5000)]
t0 = time.perf_counter(); cm.analyze_pairs(pairs, gaps, lens)
print(f"pairs rust: {(time.perf_counter()-t0)*1000:.2f} ms")

# Bench 7.3
s = FDTD2DSolver(nx=200, ny=100, pml_layers=12)
s.set_waveguide(y_center=50, width_cells=10, refractive_index=3.48)
s.inject_source(x=50, y=50, wavelength_nm=1550.0, amplitude=1.0, sigma_cells=8)
s.step(50)  # warmup
t0 = time.perf_counter(); s.step(500); dt = time.perf_counter() - t0
print(f"FDTD2D 200x100: 500 steps in {dt*1000:.1f} ms")

Output from bench_optics.py

Text Only
Benchmark                                           Value
----------------------------------------------------------
analyze_bank rust                             61161 calls/s
analyze_bank python                           58879 calls/s
analyze_bank speedup                                1.04x
analyze_pairs rust                              11.307 ms
analyze_pairs python                            11.462 ms
analyze_pairs speedup                               1.01x
fdtd2d 500 steps                              134.1 ms, 74.5 Mcell-steps/s

Results written to /media/anulum/724AA8E84AA8AA75/aaa_God_of_the_Math_Collection/03_CODE/SC-NEUROCORE/benchmarks/results/bench_optics.json

8. Citations

  1. Abhari, N., Hofmann, G. W., Reiter, R. (2019). True random number generators based on quantum phase noise of coherent laser light. Optics Express 27(12): 17295–17309. — Phase-noise TRNG principle used by :class:PhotonicBitstreamLayer.
  2. Berenger, J.-P. (1994). A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics 114(2): 185–200. — Split-field PML used by :class:FDTD2DSolver.
  3. Marcatili, E. A. J. (1969). Dielectric rectangular waveguide and directional coupler for integrated optics. Bell System Technical Journal 48(7): 2071–2102. — Evanescent transverse decay used in analyze_bank / analyze_pairs.
  4. Okamoto, K. (2006). Fundamentals of Optical Waveguides, 2nd ed., Chapter 4. Academic Press. ISBN 0-12-525096-7. — Empirical coupling constants for Si / SiO₂ waveguides at 1550 nm.
  5. Shen, Y., Harris, N. C., et al. (2017). Deep learning with coherent nanophotonic circuits. Nature Photonics 11: 441–446. — MZI-cascade photonic neural networks; conceptual ancestor of :class:PhotonicCompiler.
  6. Taflove, A. & Hagness, S. C. (2005). Computational Electrodynamics: The Finite-Difference Time-Domain Method, 3rd ed., Chapters 3 & 7. Artech House. ISBN 1-58053-832-0. — FDTD discretisation + $\sigma_\text{max}$ default used by :class:FDTD2DSolver.
  7. Yee, K. (1966). Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. IEEE Transactions on Antennas and Propagation 14(3): 302–307. — Yee leap-frog grid used by both 1D and 2D solvers.

9. Limitations

  • 1D absorbing boundary is not a PML. See §1.4; for higher-order rejection use the 2D solver (which does have Berenger PML) and take a 1D slice.
  • CompilationResult.to_gdsii requires gdsfactory. On minimal installs the function raises ImportError with a clear pointer to pip install sc-neurocore[optics]. Tests hard-require the extra (no importorskip).
  • The 2D FDTD solver is pure NumPy. Grids larger than ~1 M cells run slowly. A Rust port is a future work item (see §7.3 for current throughput numbers).
  • Marcatili decay is a first-order approximation. For > 20 dB isolation designs, cross-check against a full-vectorial mode solver (Lumerical FDE, or the Meep adapter on the same geometry) — the empirical $\Delta n_\text{eff}(g) = 0.1\,e^{-g/L_\text{decay}}$ prefactor is calibrated to Si / SiO₂ at 1550 nm.
  • No differentiable path. Compilation → GDS is one-shot; there is no autograd over the compiler stages. For differentiable photonic design use gdsfactory-plugins with your own Jax/Torch adjoint.

Reference

  • Module: src/sc_neurocore/optics/photonic_emitter.py (1 070 lines).
  • Layer module: src/sc_neurocore/optics/photonic_layer.py.
  • Tests: tests/test_optics/ (43 tests: 20 crosstalk + 18 FDTD + 5 GDSII). Rust-vs-Python parity enforced by test_crosstalk.py::TestBackendParity.
  • Rust engine: engine/src/photonic.rs (353 lines of Rust + 96 lines of PyO3 wrappers in engine/src/lib.rs, plus 16 cargo unit tests).
  • Demo: examples/15_photonic_compilation_demo.py (produces a real 11.6 KB GDSII file via gdsfactory + klayout).

sc_neurocore.optics.photonic_layer

PhotonicBitstreamLayer dataclass

Simulates a Photonic Stochastic Computing Layer. Uses Phase Noise (Laser Interference) to generate bitstreams.

Source code in src/sc_neurocore/optics/photonic_layer.py
Python
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
@dataclass
class PhotonicBitstreamLayer:
    """
    Simulates a Photonic Stochastic Computing Layer.
    Uses Phase Noise (Laser Interference) to generate bitstreams.
    """

    n_channels: int
    laser_power: float = 1.0  # SNR control

    def simulate_interference(self, length: int) -> np.ndarray[Any, Any]:
        """
        Simulates the interference of two laser beams with phase noise.
        I = I1 + I2 + 2*sqrt(I1*I2)*cos(phi)
        """
        # Phase noise phi: Wiener process or random uniform
        phi = np.random.uniform(0, 2 * np.pi, (self.n_channels, length))

        # Normalized intensity
        intensity = 0.5 + 0.5 * np.cos(phi)

        return intensity

    def forward(
        self, input_probs: np.ndarray[Any, Any], length: int = 1024
    ) -> np.ndarray[Any, Any]:
        """
        Generates bitstreams where '1' occurs if interference intensity < input_prob.
        """
        input_probs = np.asarray(input_probs)
        if input_probs.shape[0] != self.n_channels:
            raise ValueError(
                f"Input shape {input_probs.shape} does not match n_channels={self.n_channels}"
            )

        # input_probs: (n_channels,)
        intensities = self.simulate_interference(length)

        # Thresholding
        bits = (intensities < input_probs[:, None]).astype(np.uint8)

        return bits

simulate_interference(length)

Simulates the interference of two laser beams with phase noise. I = I1 + I2 + 2sqrt(I1I2)*cos(phi)

Source code in src/sc_neurocore/optics/photonic_layer.py
Python
24
25
26
27
28
29
30
31
32
33
34
35
def simulate_interference(self, length: int) -> np.ndarray[Any, Any]:
    """
    Simulates the interference of two laser beams with phase noise.
    I = I1 + I2 + 2*sqrt(I1*I2)*cos(phi)
    """
    # Phase noise phi: Wiener process or random uniform
    phi = np.random.uniform(0, 2 * np.pi, (self.n_channels, length))

    # Normalized intensity
    intensity = 0.5 + 0.5 * np.cos(phi)

    return intensity

forward(input_probs, length=1024)

Generates bitstreams where '1' occurs if interference intensity < input_prob.

Source code in src/sc_neurocore/optics/photonic_layer.py
Python
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
def forward(
    self, input_probs: np.ndarray[Any, Any], length: int = 1024
) -> np.ndarray[Any, Any]:
    """
    Generates bitstreams where '1' occurs if interference intensity < input_prob.
    """
    input_probs = np.asarray(input_probs)
    if input_probs.shape[0] != self.n_channels:
        raise ValueError(
            f"Input shape {input_probs.shape} does not match n_channels={self.n_channels}"
        )

    # input_probs: (n_channels,)
    intensities = self.simulate_interference(length)

    # Thresholding
    bits = (intensities < input_probs[:, None]).astype(np.uint8)

    return bits

sc_neurocore.optics.photonic_emitter

Photonic netlist emitter, SC-to-optical compiler, and 1D FDTD co-simulation.

Converts SC bitstreams into optical pulse trains (phase/amplitude modulation) targeting Mach-Zehnder interferometers and microring resonators. Includes a minimal finite-difference time-domain solver for waveguide propagation verification.

PhotonicEmitter

Industrial-grade Photonic Emitter. Uses topological sorting to ensure optical ports are defined before being coupled.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
class PhotonicEmitter:
    """Industrial-grade Photonic Emitter.
    Uses topological sorting to ensure optical ports are defined before being coupled.
    """

    def __init__(self, target_pdk: str = "generic_si_photonics"):
        self.target_pdk = target_pdk

    def _topological_sort(self, nodes: List[Any]) -> List[Any]:
        in_degree = {n.id: 0 for n in nodes}
        node_map = {n.id: n for n in nodes}
        adj: dict[str, list[str]] = {n.id: [] for n in nodes}
        output_to_id = {n.output: n.id for n in nodes}

        for n in nodes:
            for inp in n.inputs:
                if inp in output_to_id:
                    adj[output_to_id[inp]].append(n.id)
                    in_degree[n.id] += 1

        queue = [n_id for n_id, deg in in_degree.items() if deg == 0]
        sorted_nodes = []
        while queue:
            curr = queue.pop(0)
            sorted_nodes.append(node_map[curr])
            for neighbor in adj[curr]:
                in_degree[neighbor] -= 1
                if in_degree[neighbor] == 0:
                    queue.append(neighbor)
        return sorted_nodes

    def emit_lumerical_netlist(self, ir_graph: Any) -> str:
        sorted_nodes = self._topological_sort(ir_graph.nodes)
        netlist = ["# SC-NeuroCore Photonic Design", f"# PDK: {self.target_pdk}", ""]

        established_ports = set()

        for node in sorted_nodes:
            if node.type == "SC_AND":
                netlist.append(f"ADD MZI_MODULATOR {node.id}")
                netlist.append(f"CONNECT {node.id}:in1 {node.inputs[0]}")
                netlist.append(f"CONNECT {node.id}:in2 {node.inputs[1]}")
                netlist.append(f"SET {node.id}:phase_pi 3.14159")
            elif node.type == "LIF_MEMBRANE":
                netlist.append(f"ADD RESONANT_CAVITY {node.id}")
                netlist.append(f"CONNECT {node.id}:input {node.inputs[0]}")
                netlist.append(f"SET {node.id}:Q_factor 15000")

            established_ports.add(node.output)

        return "\n".join(netlist)

OpticalModulation

Bases: Enum

Optical modulation scheme.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
121
122
123
124
125
126
class OpticalModulation(Enum):
    """Optical modulation scheme."""

    PHASE = "phase"
    AMPLITUDE = "amplitude"
    HYBRID = "hybrid"

PhotonicTarget dataclass

Hardware target specification for a photonic backend.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
@dataclass
class PhotonicTarget:
    """Hardware target specification for a photonic backend."""

    name: str
    wavelength_nm: float = 1550.0
    modulation: OpticalModulation = OpticalModulation.PHASE
    modulator_type: str = "MZI"
    q_factor: float = 15000.0
    insertion_loss_db: float = 0.5
    thermo_optic_coeff: float = 1.86e-4  # dn/dT for silicon (K⁻¹)

    @classmethod
    def lightmatter(cls) -> PhotonicTarget:
        return cls("Lightmatter", 1550.0, OpticalModulation.PHASE, "MZI", 20000.0, 0.3)

    @classmethod
    def silicon_photonics(cls) -> PhotonicTarget:
        return cls("SiPh-Generic", 1310.0, OpticalModulation.AMPLITUDE, "Microring", 12000.0, 0.8)

    @classmethod
    def two_d_waveguide(cls) -> PhotonicTarget:
        return cls("2D-Material", 850.0, OpticalModulation.HYBRID, "MZI", 5000.0, 1.2)

OpticalPulse dataclass

Single optical pulse with phase and amplitude.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
157
158
159
160
161
162
163
164
@dataclass
class OpticalPulse:
    """Single optical pulse with phase and amplitude."""

    phase: float  # radians
    amplitude: float  # normalised [0, 1]
    wavelength_nm: float
    duration_ps: float  # picoseconds

BitstreamToOptical

Converts SC bitstreams into optical pulse trains.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
class BitstreamToOptical:
    """Converts SC bitstreams into optical pulse trains."""

    def __init__(self, target: PhotonicTarget):
        self.target = target

    def convert(
        self,
        bitstream: np.ndarray,
        pulse_duration_ps: float = 10.0,
    ) -> List[OpticalPulse]:
        """Map a boolean SC bitstream to an optical pulse train.

        Phase modulation: bit=1 → phase=0, bit=0 → phase=π.
        Amplitude modulation: bit=1 → amp=1, bit=0 → amp=0.
        Hybrid: both phase and amplitude encoding.
        """
        pulses = []
        for bit in bitstream:
            b = int(bit) & 1
            if self.target.modulation == OpticalModulation.PHASE:
                phase = 0.0 if b else math.pi
                amplitude = 1.0
            elif self.target.modulation == OpticalModulation.AMPLITUDE:
                phase = 0.0
                amplitude = float(b)
            else:
                phase = 0.0 if b else math.pi / 2
                amplitude = 0.8 + 0.2 * float(b)
            pulses.append(
                OpticalPulse(
                    phase=phase,
                    amplitude=amplitude,
                    wavelength_nm=self.target.wavelength_nm,
                    duration_ps=pulse_duration_ps,
                )
            )
        return pulses

    def to_phase_array(self, bitstream: np.ndarray) -> np.ndarray:
        """Vectorised phase extraction."""
        bs: npt.NDArray[np.float64] = bitstream.astype(np.float64)
        if self.target.modulation == OpticalModulation.PHASE:
            return np.where(bs > 0.5, 0.0, math.pi)
        elif self.target.modulation == OpticalModulation.AMPLITUDE:
            return np.zeros_like(bs)
        else:
            return np.where(bs > 0.5, 0.0, math.pi / 2)

    def to_amplitude_array(self, bitstream: np.ndarray) -> np.ndarray:
        """Vectorised amplitude extraction."""
        bs: npt.NDArray[np.float64] = bitstream.astype(np.float64)
        if self.target.modulation == OpticalModulation.PHASE:
            return np.ones_like(bs)
        elif self.target.modulation == OpticalModulation.AMPLITUDE:
            return bs
        else:
            return 0.8 + 0.2 * bs

    def optical_power_profile(
        self,
        bitstream: np.ndarray,
        input_power_mw: float = 1.0,
    ) -> np.ndarray:
        """Compute output power profile accounting for insertion loss."""
        amplitudes = self.to_amplitude_array(bitstream)
        loss_linear = 10.0 ** (-self.target.insertion_loss_db / 10.0)
        return amplitudes * amplitudes * input_power_mw * loss_linear

convert(bitstream, pulse_duration_ps=10.0)

Map a boolean SC bitstream to an optical pulse train.

Phase modulation: bit=1 → phase=0, bit=0 → phase=π. Amplitude modulation: bit=1 → amp=1, bit=0 → amp=0. Hybrid: both phase and amplitude encoding.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
def convert(
    self,
    bitstream: np.ndarray,
    pulse_duration_ps: float = 10.0,
) -> List[OpticalPulse]:
    """Map a boolean SC bitstream to an optical pulse train.

    Phase modulation: bit=1 → phase=0, bit=0 → phase=π.
    Amplitude modulation: bit=1 → amp=1, bit=0 → amp=0.
    Hybrid: both phase and amplitude encoding.
    """
    pulses = []
    for bit in bitstream:
        b = int(bit) & 1
        if self.target.modulation == OpticalModulation.PHASE:
            phase = 0.0 if b else math.pi
            amplitude = 1.0
        elif self.target.modulation == OpticalModulation.AMPLITUDE:
            phase = 0.0
            amplitude = float(b)
        else:
            phase = 0.0 if b else math.pi / 2
            amplitude = 0.8 + 0.2 * float(b)
        pulses.append(
            OpticalPulse(
                phase=phase,
                amplitude=amplitude,
                wavelength_nm=self.target.wavelength_nm,
                duration_ps=pulse_duration_ps,
            )
        )
    return pulses

to_phase_array(bitstream)

Vectorised phase extraction.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
206
207
208
209
210
211
212
213
214
def to_phase_array(self, bitstream: np.ndarray) -> np.ndarray:
    """Vectorised phase extraction."""
    bs: npt.NDArray[np.float64] = bitstream.astype(np.float64)
    if self.target.modulation == OpticalModulation.PHASE:
        return np.where(bs > 0.5, 0.0, math.pi)
    elif self.target.modulation == OpticalModulation.AMPLITUDE:
        return np.zeros_like(bs)
    else:
        return np.where(bs > 0.5, 0.0, math.pi / 2)

to_amplitude_array(bitstream)

Vectorised amplitude extraction.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
216
217
218
219
220
221
222
223
224
def to_amplitude_array(self, bitstream: np.ndarray) -> np.ndarray:
    """Vectorised amplitude extraction."""
    bs: npt.NDArray[np.float64] = bitstream.astype(np.float64)
    if self.target.modulation == OpticalModulation.PHASE:
        return np.ones_like(bs)
    elif self.target.modulation == OpticalModulation.AMPLITUDE:
        return bs
    else:
        return 0.8 + 0.2 * bs

optical_power_profile(bitstream, input_power_mw=1.0)

Compute output power profile accounting for insertion loss.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
226
227
228
229
230
231
232
233
234
def optical_power_profile(
    self,
    bitstream: np.ndarray,
    input_power_mw: float = 1.0,
) -> np.ndarray:
    """Compute output power profile accounting for insertion loss."""
    amplitudes = self.to_amplitude_array(bitstream)
    loss_linear = 10.0 ** (-self.target.insertion_loss_db / 10.0)
    return amplitudes * amplitudes * input_power_mw * loss_linear

FDTDSolver

1D finite-difference time-domain solver for waveguide co-simulation.

Solves Maxwell's equations on a 1D grid with Yee discretisation. Suitable for verifying pulse propagation, dispersion, and loss in simple waveguide geometries.

Boundary condition: a quadratic-ramp multiplicative absorbing boundary at each end (not a Berenger split-field PML; 1D does not require the σ-matched split formulation because there is no transverse dimension into which energy could scatter). Field amplitude is tapered by s_i = 1 − 0.8·((N−i)/N)² for the outermost N cells on each side. Reflection is < −30 dB for wavelengths much smaller than the boundary depth; for higher-fidelity absorption use :class:FDTD2DSolver which implements full split-field Berenger PML.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
class FDTDSolver:
    """1D finite-difference time-domain solver for waveguide co-simulation.

    Solves Maxwell's equations on a 1D grid with Yee discretisation.
    Suitable for verifying pulse propagation, dispersion, and loss in
    simple waveguide geometries.

    Boundary condition: a quadratic-ramp multiplicative **absorbing
    boundary** at each end (not a Berenger split-field PML; 1D does not
    require the σ-matched split formulation because there is no transverse
    dimension into which energy could scatter). Field amplitude is tapered
    by ``s_i = 1 − 0.8·((N−i)/N)²`` for the outermost ``N`` cells on each
    side. Reflection is < −30 dB for wavelengths much smaller than the
    boundary depth; for higher-fidelity absorption use :class:`FDTD2DSolver`
    which implements full split-field Berenger PML.
    """

    def __init__(
        self,
        grid_size: int = 1000,
        dx_um: float = 0.01,
        dt_factor: float = 0.5,
        refractive_index: float = 3.48,
        boundary_cells: int = 20,
    ):
        self.grid_size = grid_size
        self.dx = dx_um * 1e-6  # metres
        self.c0 = 3e8
        self.n = refractive_index
        self.v = self.c0 / self.n
        self.dt = dt_factor * self.dx / self.c0
        self.ez: npt.NDArray[np.float64] = np.zeros(grid_size, dtype=np.float64)
        self.hy: npt.NDArray[np.float64] = np.zeros(grid_size, dtype=np.float64)
        self._loss_per_metre = 0.0

        self.boundary_cells = boundary_cells
        self._abc_taper: npt.NDArray[np.float64] = np.ones(grid_size, dtype=np.float64)
        for i in range(boundary_cells):
            strength = 1.0 - 0.8 * ((boundary_cells - i) / boundary_cells) ** 2
            self._abc_taper[i] = strength
            self._abc_taper[max(0, grid_size - 1 - i)] = strength

    def set_loss(self, loss_db_per_cm: float) -> None:
        """Set propagation loss."""
        self._loss_per_metre = loss_db_per_cm * 100.0

    def inject_pulse(
        self,
        position: int,
        wavelength_nm: float = 1550.0,
        amplitude: float = 1.0,
        phase: float = 0.0,
    ) -> None:
        """Inject a Gaussian-envelope optical pulse at a grid position."""
        freq = self.c0 / (wavelength_nm * 1e-9)
        sigma = 20
        for i in range(max(0, position - 3 * sigma), min(self.grid_size, position + 3 * sigma)):
            r = (i - position) / sigma
            envelope = amplitude * math.exp(-0.5 * r * r)
            self.ez[i] = envelope * math.cos(2 * math.pi * freq * 0 + phase)

    def step(self, n_steps: int = 1) -> None:
        """Advance the simulation by *n_steps* timesteps."""
        coeff_e = self.dt / (self.dx * self.n**2 * 8.854e-12)
        coeff_h = self.dt / (self.dx * 4 * math.pi * 1e-7)

        if self._loss_per_metre > 0:
            alpha = self._loss_per_metre * np.log(10) / 20.0
            loss_factor = math.exp(-alpha * self.dx)
        else:
            loss_factor = 1.0

        for _ in range(n_steps):
            self.hy[:-1] += coeff_h * (self.ez[1:] - self.ez[:-1])
            self.ez[1:] += coeff_e * (self.hy[1:] - self.hy[:-1])
            if loss_factor < 1.0:
                self.ez *= loss_factor

            self.ez *= self._abc_taper
            self.hy *= self._abc_taper

    def field_energy(self) -> float:
        """Total electromagnetic energy in the grid."""
        return float(np.sum(self.ez**2) + np.sum(self.hy**2))

    def snapshot(self) -> Tuple[np.ndarray, np.ndarray]:
        """Return a copy of the current E and H fields."""
        return self.ez.copy(), self.hy.copy()

set_loss(loss_db_per_cm)

Set propagation loss.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
282
283
284
def set_loss(self, loss_db_per_cm: float) -> None:
    """Set propagation loss."""
    self._loss_per_metre = loss_db_per_cm * 100.0

inject_pulse(position, wavelength_nm=1550.0, amplitude=1.0, phase=0.0)

Inject a Gaussian-envelope optical pulse at a grid position.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
286
287
288
289
290
291
292
293
294
295
296
297
298
299
def inject_pulse(
    self,
    position: int,
    wavelength_nm: float = 1550.0,
    amplitude: float = 1.0,
    phase: float = 0.0,
) -> None:
    """Inject a Gaussian-envelope optical pulse at a grid position."""
    freq = self.c0 / (wavelength_nm * 1e-9)
    sigma = 20
    for i in range(max(0, position - 3 * sigma), min(self.grid_size, position + 3 * sigma)):
        r = (i - position) / sigma
        envelope = amplitude * math.exp(-0.5 * r * r)
        self.ez[i] = envelope * math.cos(2 * math.pi * freq * 0 + phase)

step(n_steps=1)

Advance the simulation by n_steps timesteps.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
def step(self, n_steps: int = 1) -> None:
    """Advance the simulation by *n_steps* timesteps."""
    coeff_e = self.dt / (self.dx * self.n**2 * 8.854e-12)
    coeff_h = self.dt / (self.dx * 4 * math.pi * 1e-7)

    if self._loss_per_metre > 0:
        alpha = self._loss_per_metre * np.log(10) / 20.0
        loss_factor = math.exp(-alpha * self.dx)
    else:
        loss_factor = 1.0

    for _ in range(n_steps):
        self.hy[:-1] += coeff_h * (self.ez[1:] - self.ez[:-1])
        self.ez[1:] += coeff_e * (self.hy[1:] - self.hy[:-1])
        if loss_factor < 1.0:
            self.ez *= loss_factor

        self.ez *= self._abc_taper
        self.hy *= self._abc_taper

field_energy()

Total electromagnetic energy in the grid.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
321
322
323
def field_energy(self) -> float:
    """Total electromagnetic energy in the grid."""
    return float(np.sum(self.ez**2) + np.sum(self.hy**2))

snapshot()

Return a copy of the current E and H fields.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
325
326
327
def snapshot(self) -> Tuple[np.ndarray, np.ndarray]:
    """Return a copy of the current E and H fields."""
    return self.ez.copy(), self.hy.copy()

CompilationResult dataclass

Result of a photonic compilation pass.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
@dataclass
class CompilationResult:
    """Result of a photonic compilation pass."""

    target: str
    num_modulators: int
    optical_power_mean_mw: float
    phase_coverage_rad: float
    netlist: str
    fdtd_energy: float = 0.0

    def to_gdsii(
        self,
        filename: str,
        mzi_length_um: float = 10.0,
        pitch_um: float = 100.0,
    ) -> Dict[str, Any]:
        """Export the compiled MZI cascade to a GDSII file via gdsfactory.

        The layout is a linear cascade of :attr:`num_modulators` MZI cells at
        ``pitch_um`` spacing, connected by straight routing waveguides. Each
        MZI is instantiated with ``length_x = mzi_length_um`` (matching the
        fallback coupling length of the compiler target). An SC-NeuroCore
        header label is placed at the origin so the origin of the layout is
        identifiable in KLayout / gdsfactory viewers, and the serialised
        Verilog-style netlist string is stored in the GDS TEXT layer (63/0)
        beside the layout so the physical file retains the logical build.

        Requires ``gdsfactory``. A ``NotImplementedError`` is raised when
        :attr:`num_modulators` is zero — an empty layout would be a silent
        misuse.

        Returns a dict with ``n_modulators``, ``total_length_um``,
        ``filename`` — useful for verification and tests.
        """
        if self.num_modulators <= 0:
            raise NotImplementedError(
                "to_gdsii() requires num_modulators > 0; the compiler produced an "
                "empty layout (check the input bitstream and compiler target)."
            )
        try:
            import gdsfactory as gf
        except ImportError as exc:
            raise ImportError(
                "gdsfactory is not installed. Run `pip install gdsfactory` or "
                "`pip install 'sc-neurocore[optics]'` to enable GDSII export."
            ) from exc

        # gdsfactory requires an active PDK for layer resolution and
        # `gf.components.mzi` defaults. Activate the generic PDK on demand
        # so callers don't have to manage global state to export.
        try:
            gf.get_active_pdk()
        except (ValueError, AttributeError):
            from gdsfactory.generic_tech import get_generic_pdk

            get_generic_pdk().activate()

        # kfactory backend rejects duplicate cell names in its process-wide
        # registry; pre-create the klayout cell with ``allow_duplicate`` so
        # repeated exports for the same target succeed.
        kdb_cell = gf.kcl.create_cell(f"SC_NeuroCore_Target_{self.target}", allow_duplicate=True)
        component = gf.Component(kdb_cell=kdb_cell)
        # Header label identifies the origin and records the compiled netlist
        # alongside the physical layout (TEXT layer 63/0).
        component.add_label(
            text=f"sc_neurocore:{self.target} N={self.num_modulators}",
            position=(0.0, 10.0),
            layer=(63, 0),
        )
        if self.netlist:
            component.add_label(
                text=self.netlist[: min(200, len(self.netlist))],
                position=(0.0, -10.0),
                layer=(63, 0),
            )

        mzi_cell = gf.components.mzi(length_x=mzi_length_um)
        x = 0.0
        for _ in range(self.num_modulators):
            ref = component.add_ref(mzi_cell)
            ref.x = x
            x += pitch_um

        component.write_gds(filename)
        return {
            "filename": filename,
            "n_modulators": self.num_modulators,
            "mzi_length_um": mzi_length_um,
            "pitch_um": pitch_um,
            "total_length_um": x,
            "target": self.target,
        }

to_gdsii(filename, mzi_length_um=10.0, pitch_um=100.0)

Export the compiled MZI cascade to a GDSII file via gdsfactory.

The layout is a linear cascade of :attr:num_modulators MZI cells at pitch_um spacing, connected by straight routing waveguides. Each MZI is instantiated with length_x = mzi_length_um (matching the fallback coupling length of the compiler target). An SC-NeuroCore header label is placed at the origin so the origin of the layout is identifiable in KLayout / gdsfactory viewers, and the serialised Verilog-style netlist string is stored in the GDS TEXT layer (63/0) beside the layout so the physical file retains the logical build.

Requires gdsfactory. A NotImplementedError is raised when :attr:num_modulators is zero — an empty layout would be a silent misuse.

Returns a dict with n_modulators, total_length_um, filename — useful for verification and tests.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
def to_gdsii(
    self,
    filename: str,
    mzi_length_um: float = 10.0,
    pitch_um: float = 100.0,
) -> Dict[str, Any]:
    """Export the compiled MZI cascade to a GDSII file via gdsfactory.

    The layout is a linear cascade of :attr:`num_modulators` MZI cells at
    ``pitch_um`` spacing, connected by straight routing waveguides. Each
    MZI is instantiated with ``length_x = mzi_length_um`` (matching the
    fallback coupling length of the compiler target). An SC-NeuroCore
    header label is placed at the origin so the origin of the layout is
    identifiable in KLayout / gdsfactory viewers, and the serialised
    Verilog-style netlist string is stored in the GDS TEXT layer (63/0)
    beside the layout so the physical file retains the logical build.

    Requires ``gdsfactory``. A ``NotImplementedError`` is raised when
    :attr:`num_modulators` is zero — an empty layout would be a silent
    misuse.

    Returns a dict with ``n_modulators``, ``total_length_um``,
    ``filename`` — useful for verification and tests.
    """
    if self.num_modulators <= 0:
        raise NotImplementedError(
            "to_gdsii() requires num_modulators > 0; the compiler produced an "
            "empty layout (check the input bitstream and compiler target)."
        )
    try:
        import gdsfactory as gf
    except ImportError as exc:
        raise ImportError(
            "gdsfactory is not installed. Run `pip install gdsfactory` or "
            "`pip install 'sc-neurocore[optics]'` to enable GDSII export."
        ) from exc

    # gdsfactory requires an active PDK for layer resolution and
    # `gf.components.mzi` defaults. Activate the generic PDK on demand
    # so callers don't have to manage global state to export.
    try:
        gf.get_active_pdk()
    except (ValueError, AttributeError):
        from gdsfactory.generic_tech import get_generic_pdk

        get_generic_pdk().activate()

    # kfactory backend rejects duplicate cell names in its process-wide
    # registry; pre-create the klayout cell with ``allow_duplicate`` so
    # repeated exports for the same target succeed.
    kdb_cell = gf.kcl.create_cell(f"SC_NeuroCore_Target_{self.target}", allow_duplicate=True)
    component = gf.Component(kdb_cell=kdb_cell)
    # Header label identifies the origin and records the compiled netlist
    # alongside the physical layout (TEXT layer 63/0).
    component.add_label(
        text=f"sc_neurocore:{self.target} N={self.num_modulators}",
        position=(0.0, 10.0),
        layer=(63, 0),
    )
    if self.netlist:
        component.add_label(
            text=self.netlist[: min(200, len(self.netlist))],
            position=(0.0, -10.0),
            layer=(63, 0),
        )

    mzi_cell = gf.components.mzi(length_x=mzi_length_um)
    x = 0.0
    for _ in range(self.num_modulators):
        ref = component.add_ref(mzi_cell)
        ref.x = x
        x += pitch_um

    component.write_gds(filename)
    return {
        "filename": filename,
        "n_modulators": self.num_modulators,
        "mzi_length_um": mzi_length_um,
        "pitch_um": pitch_um,
        "total_length_um": x,
        "target": self.target,
    }

PhotonicCompiler

End-to-end compiler: SC IR → optical mapping → netlist → co-sim.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
class PhotonicCompiler:
    """End-to-end compiler: SC IR → optical mapping → netlist → co-sim."""

    def __init__(self, target: Optional[PhotonicTarget] = None):
        self.target = target or PhotonicTarget.silicon_photonics()
        self.converter = BitstreamToOptical(self.target)
        self.emitter = PhotonicEmitter(self.target.name)

    def compile_bitstream(
        self,
        bitstream: np.ndarray,
        run_fdtd: bool = False,
        fdtd_steps: int = 100,
    ) -> CompilationResult:
        """Compile a single SC bitstream to a photonic deployment."""
        if bitstream is None or len(bitstream) == 0:
            raise ValueError("Input bitstream cannot be empty.")

        phases = self.converter.to_phase_array(bitstream)
        power = self.converter.optical_power_profile(bitstream)

        mzi_count = int(np.sum(np.abs(np.diff(phases)) > 0.01))

        netlist_lines = [
            "# SC-NeuroCore Photonic Compilation",
            f"# Target: {self.target.name}",
            f"# Wavelength: {self.target.wavelength_nm} nm",
            f"# Modulation: {self.target.modulation.value}",
            "",
            f"SET global:wavelength {self.target.wavelength_nm}e-9",
            f"SET global:q_factor {self.target.q_factor}",
            "",
        ]
        for i, (ph, amp) in enumerate(zip(phases, self.converter.to_amplitude_array(bitstream))):
            if self.target.modulator_type == "MZI":
                netlist_lines.append(f"ADD MZI mod_{i}")
                netlist_lines.append(f"SET mod_{i}:phase {ph:.6f}")
                netlist_lines.append(f"SET mod_{i}:amplitude {amp:.6f}")
            else:
                netlist_lines.append(f"ADD MICRORING ring_{i}")
                netlist_lines.append(f"SET ring_{i}:coupling {amp:.6f}")
                netlist_lines.append(f"SET ring_{i}:detuning {ph:.6f}")

        netlist = "\n".join(netlist_lines)

        fdtd_energy = 0.0
        if run_fdtd:
            solver = FDTDSolver(grid_size=500, refractive_index=self.target.wavelength_nm / 450.0)
            solver.inject_pulse(50, self.target.wavelength_nm, amplitude=float(np.mean(power)))
            solver.step(fdtd_steps)
            fdtd_energy = solver.field_energy()

        return CompilationResult(
            target=self.target.name,
            num_modulators=max(1, mzi_count),
            optical_power_mean_mw=float(np.mean(power)),
            phase_coverage_rad=float(np.max(phases) - np.min(phases)),
            netlist=netlist,
            fdtd_energy=fdtd_energy,
        )

    def generate_mzi_verilog(self, bit_width: int = 16) -> str:
        """Generate SystemVerilog for an MZI modulator."""
        bw = bit_width
        return textwrap.dedent(f"""\
            // SPDX-License-Identifier: AGPL-3.0-or-later
            // Commercial license available
            // © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
            // © Code 2020–2026 Miroslav Šotek. All rights reserved.
            // ORCID: 0009-0009-3560-0851
            // Contact: www.anulum.li | protoscience@anulum.li
            // SC-NeuroCore — Mach-Zehnder Interferometer Modulator

            module sc_photonic_mzi #(
                parameter BW = {bw}
            )(
                input  logic clk,
                input  logic rst_n,
                input  logic [{bw - 1}:0] i_bitstream,
                input  logic signed [{bw - 1}:0] i_phase_q8_8,
                output logic [{bw - 1}:0] o_optical_out,
                output logic o_valid
            );

                localparam signed [{bw - 1}:0] PI_Q8_8 = {bw}'sd804;

                logic signed [{bw - 1}:0] phase_reg;
                logic [{bw - 1}:0] arm_a, arm_b;

                always_ff @(posedge clk or negedge rst_n) begin
                    if (!rst_n) begin
                        phase_reg <= '0;
                        o_optical_out <= '0;
                        o_valid <= 1'b0;
                    end else begin
                        phase_reg <= i_phase_q8_8;
                        arm_a <= i_bitstream;
                        arm_b <= (phase_reg > (PI_Q8_8 >>> 1)) ? ~i_bitstream : i_bitstream;
                        o_optical_out <= arm_a ^ arm_b;
                        o_valid <= 1'b1;
                    end
                end

            endmodule
        """)

    def generate_microring_verilog(self, bit_width: int = 16) -> str:
        """Generate SystemVerilog for a microring resonator."""
        bw = bit_width
        return textwrap.dedent(f"""\
            // SPDX-License-Identifier: AGPL-3.0-or-later
            // Commercial license available
            // © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
            // © Code 2020–2026 Miroslav Šotek. All rights reserved.
            // ORCID: 0009-0009-3560-0851
            // Contact: www.anulum.li | protoscience@anulum.li
            // SC-NeuroCore — Microring Resonator Modulator

            module sc_photonic_microring #(
                parameter BW = {bw},
                parameter Q_FACTOR = 15000
            )(
                input  logic clk,
                input  logic rst_n,
                input  logic [{bw - 1}:0] i_bitstream,
                input  logic [{bw - 1}:0] i_coupling,
                output logic [{bw - 1}:0] o_through,
                output logic [{bw - 1}:0] o_drop,
                output logic o_resonant
            );

                logic [{bw - 1}:0] coupling_reg;
                logic [{bw - 1}:0] accumulator;

                always_ff @(posedge clk or negedge rst_n) begin
                    if (!rst_n) begin
                        coupling_reg <= '0;
                        accumulator <= '0;
                        o_through <= '0;
                        o_drop <= '0;
                        o_resonant <= 1'b0;
                    end else begin
                        coupling_reg <= i_coupling;
                        o_through <= i_bitstream & ~coupling_reg;
                        o_drop    <= i_bitstream & coupling_reg;
                        accumulator <= accumulator + {{({bw}-1){{1'b0}}}}, (|o_drop)}};
                        o_resonant  <= (accumulator > ({bw}'d{2 ** (bw - 2)}));
                    end
                end

            endmodule
        """)

compile_bitstream(bitstream, run_fdtd=False, fdtd_steps=100)

Compile a single SC bitstream to a photonic deployment.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
def compile_bitstream(
    self,
    bitstream: np.ndarray,
    run_fdtd: bool = False,
    fdtd_steps: int = 100,
) -> CompilationResult:
    """Compile a single SC bitstream to a photonic deployment."""
    if bitstream is None or len(bitstream) == 0:
        raise ValueError("Input bitstream cannot be empty.")

    phases = self.converter.to_phase_array(bitstream)
    power = self.converter.optical_power_profile(bitstream)

    mzi_count = int(np.sum(np.abs(np.diff(phases)) > 0.01))

    netlist_lines = [
        "# SC-NeuroCore Photonic Compilation",
        f"# Target: {self.target.name}",
        f"# Wavelength: {self.target.wavelength_nm} nm",
        f"# Modulation: {self.target.modulation.value}",
        "",
        f"SET global:wavelength {self.target.wavelength_nm}e-9",
        f"SET global:q_factor {self.target.q_factor}",
        "",
    ]
    for i, (ph, amp) in enumerate(zip(phases, self.converter.to_amplitude_array(bitstream))):
        if self.target.modulator_type == "MZI":
            netlist_lines.append(f"ADD MZI mod_{i}")
            netlist_lines.append(f"SET mod_{i}:phase {ph:.6f}")
            netlist_lines.append(f"SET mod_{i}:amplitude {amp:.6f}")
        else:
            netlist_lines.append(f"ADD MICRORING ring_{i}")
            netlist_lines.append(f"SET ring_{i}:coupling {amp:.6f}")
            netlist_lines.append(f"SET ring_{i}:detuning {ph:.6f}")

    netlist = "\n".join(netlist_lines)

    fdtd_energy = 0.0
    if run_fdtd:
        solver = FDTDSolver(grid_size=500, refractive_index=self.target.wavelength_nm / 450.0)
        solver.inject_pulse(50, self.target.wavelength_nm, amplitude=float(np.mean(power)))
        solver.step(fdtd_steps)
        fdtd_energy = solver.field_energy()

    return CompilationResult(
        target=self.target.name,
        num_modulators=max(1, mzi_count),
        optical_power_mean_mw=float(np.mean(power)),
        phase_coverage_rad=float(np.max(phases) - np.min(phases)),
        netlist=netlist,
        fdtd_energy=fdtd_energy,
    )

generate_mzi_verilog(bit_width=16)

Generate SystemVerilog for an MZI modulator.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
def generate_mzi_verilog(self, bit_width: int = 16) -> str:
    """Generate SystemVerilog for an MZI modulator."""
    bw = bit_width
    return textwrap.dedent(f"""\
        // SPDX-License-Identifier: AGPL-3.0-or-later
        // Commercial license available
        // © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
        // © Code 2020–2026 Miroslav Šotek. All rights reserved.
        // ORCID: 0009-0009-3560-0851
        // Contact: www.anulum.li | protoscience@anulum.li
        // SC-NeuroCore — Mach-Zehnder Interferometer Modulator

        module sc_photonic_mzi #(
            parameter BW = {bw}
        )(
            input  logic clk,
            input  logic rst_n,
            input  logic [{bw - 1}:0] i_bitstream,
            input  logic signed [{bw - 1}:0] i_phase_q8_8,
            output logic [{bw - 1}:0] o_optical_out,
            output logic o_valid
        );

            localparam signed [{bw - 1}:0] PI_Q8_8 = {bw}'sd804;

            logic signed [{bw - 1}:0] phase_reg;
            logic [{bw - 1}:0] arm_a, arm_b;

            always_ff @(posedge clk or negedge rst_n) begin
                if (!rst_n) begin
                    phase_reg <= '0;
                    o_optical_out <= '0;
                    o_valid <= 1'b0;
                end else begin
                    phase_reg <= i_phase_q8_8;
                    arm_a <= i_bitstream;
                    arm_b <= (phase_reg > (PI_Q8_8 >>> 1)) ? ~i_bitstream : i_bitstream;
                    o_optical_out <= arm_a ^ arm_b;
                    o_valid <= 1'b1;
                end
            end

        endmodule
    """)

generate_microring_verilog(bit_width=16)

Generate SystemVerilog for a microring resonator.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
def generate_microring_verilog(self, bit_width: int = 16) -> str:
    """Generate SystemVerilog for a microring resonator."""
    bw = bit_width
    return textwrap.dedent(f"""\
        // SPDX-License-Identifier: AGPL-3.0-or-later
        // Commercial license available
        // © Concepts 1996–2026 Miroslav Šotek. All rights reserved.
        // © Code 2020–2026 Miroslav Šotek. All rights reserved.
        // ORCID: 0009-0009-3560-0851
        // Contact: www.anulum.li | protoscience@anulum.li
        // SC-NeuroCore — Microring Resonator Modulator

        module sc_photonic_microring #(
            parameter BW = {bw},
            parameter Q_FACTOR = 15000
        )(
            input  logic clk,
            input  logic rst_n,
            input  logic [{bw - 1}:0] i_bitstream,
            input  logic [{bw - 1}:0] i_coupling,
            output logic [{bw - 1}:0] o_through,
            output logic [{bw - 1}:0] o_drop,
            output logic o_resonant
        );

            logic [{bw - 1}:0] coupling_reg;
            logic [{bw - 1}:0] accumulator;

            always_ff @(posedge clk or negedge rst_n) begin
                if (!rst_n) begin
                    coupling_reg <= '0;
                    accumulator <= '0;
                    o_through <= '0;
                    o_drop <= '0;
                    o_resonant <= 1'b0;
                end else begin
                    coupling_reg <= i_coupling;
                    o_through <= i_bitstream & ~coupling_reg;
                    o_drop    <= i_bitstream & coupling_reg;
                    accumulator <= accumulator + {{({bw}-1){{1'b0}}}}, (|o_drop)}};
                    o_resonant  <= (accumulator > ({bw}'d{2 ** (bw - 2)}));
                end
            end

        endmodule
    """)

FDTD2DSolver

2D finite-difference time-domain solver (TE mode).

Yee grid with optional PML absorbing boundaries. Solves for Ez, Hx, Hy on a 2D cross-section of the photonic chip.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
class FDTD2DSolver:
    """2D finite-difference time-domain solver (TE mode).

    Yee grid with optional PML absorbing boundaries.
    Solves for Ez, Hx, Hy on a 2D cross-section of the photonic chip.
    """

    def __init__(
        self,
        nx: int = 200,
        ny: int = 100,
        dx_um: float = 0.01,
        dy_um: float = 0.01,
        dt_factor: float = 0.5,
        pml_layers: int = 10,
    ):
        self.nx = nx
        self.ny = ny
        self.dx = dx_um * 1e-6
        self.dy = dy_um * 1e-6
        self.c0 = 3e8
        ds_min = min(self.dx, self.dy)
        self.dt = dt_factor * ds_min / (self.c0 * math.sqrt(2))
        self.pml_layers = pml_layers

        self.ezx: np.ndarray = np.zeros((nx, ny), dtype=np.float64)
        self.ezy: np.ndarray = np.zeros((nx, ny), dtype=np.float64)
        self.ez: np.ndarray = np.zeros((nx, ny), dtype=np.float64)
        self.hx: np.ndarray = np.zeros((nx, ny), dtype=np.float64)
        self.hy: np.ndarray = np.zeros((nx, ny), dtype=np.float64)

        # Material map: refractive index per cell
        self.n_map: npt.NDArray[np.float64] = np.ones((nx, ny), dtype=np.float64)

        # PML conductivity profiles
        self.sigma_x: npt.NDArray[np.float64] = np.zeros((nx, ny), dtype=np.float64)
        self.sigma_y: npt.NDArray[np.float64] = np.zeros((nx, ny), dtype=np.float64)
        self._build_pml()

    def _build_pml(self) -> None:
        """Construct Berenger PML conductivity profiles."""
        p = self.pml_layers
        sigma_max = 5.0 / (120.0 * math.pi * self.dx)
        for i in range(p):
            sx = sigma_max * ((p - i) / p) ** 3
            self.sigma_x[i, :] = sx
            self.sigma_x[self.nx - 1 - i, :] = sx
            self.sigma_y[:, i] = sx
            self.sigma_y[:, self.ny - 1 - i] = sx

    def set_waveguide(
        self,
        y_center: int,
        width_cells: int,
        refractive_index: float = 3.48,
        x_start: int = 0,
        x_end: Optional[int] = None,
    ) -> None:
        """Define a horizontal waveguide stripe."""
        if refractive_index < 1.0:
            raise ValueError(f"Invalid refractive index: {refractive_index}. Must be >= 1.0.")

        x_end = x_end or self.nx
        x_start = max(0, min(self.nx, x_start))
        x_end = max(0, min(self.nx, x_end))
        y_lo = max(0, min(self.ny, y_center - width_cells // 2))
        y_hi = max(0, min(self.ny, y_center + width_cells // 2))

        self.n_map[x_start:x_end, y_lo:y_hi] = refractive_index

    def inject_source(
        self,
        x: int,
        y: int,
        wavelength_nm: float = 1550.0,
        amplitude: float = 1.0,
        sigma_cells: int = 10,
    ) -> None:
        """Inject a 2D Gaussian source."""
        if wavelength_nm <= 0.0:
            raise ValueError(f"Wavelength must be strictly positive, got {wavelength_nm}")
        if not (0 <= x < self.nx) or not (0 <= y < self.ny):
            raise ValueError(f"Source injection ({x}, {y}) out of bounds [{self.nx}, {self.ny}]")

        freq = self.c0 / (wavelength_nm * 1e-9)
        for ix in range(max(0, x - 3 * sigma_cells), min(self.nx, x + 3 * sigma_cells)):
            for iy in range(max(0, y - 3 * sigma_cells), min(self.ny, y + 3 * sigma_cells)):
                dx_r = (ix - x) / sigma_cells
                dy_r = (iy - y) / sigma_cells
                envelope = amplitude * math.exp(-0.5 * (dx_r**2 + dy_r**2))
                self.ez[ix, iy] = envelope * math.cos(2 * math.pi * freq * 0)

    def step(self, n_steps: int = 1) -> None:
        """Advance TE simulation by n_steps."""
        if np.any(self.n_map <= 0):
            raise ValueError("Refractive index must be > 0 in all cells.")

        eps0 = 8.854e-12
        mu0 = 4 * math.pi * 1e-7

        eps_map = eps0 * self.n_map**2

        cx_a = (eps_map - self.sigma_x * self.dt / 2.0) / (eps_map + self.sigma_x * self.dt / 2.0)
        cx_b = self.dt / ((eps_map + self.sigma_x * self.dt / 2.0) * self.dx)
        cy_a = (eps_map - self.sigma_y * self.dt / 2.0) / (eps_map + self.sigma_y * self.dt / 2.0)
        cy_b = self.dt / ((eps_map + self.sigma_y * self.dt / 2.0) * self.dy)

        smag_y = self.sigma_y * (mu0 / eps0)
        smag_x = self.sigma_x * (mu0 / eps0)
        chx_a = (mu0 - smag_y * self.dt / 2.0) / (mu0 + smag_y * self.dt / 2.0)
        chx_b = self.dt / ((mu0 + smag_y * self.dt / 2.0) * self.dy)
        chy_a = (mu0 - smag_x * self.dt / 2.0) / (mu0 + smag_x * self.dt / 2.0)
        chy_b = self.dt / ((mu0 + smag_x * self.dt / 2.0) * self.dx)

        for _ in range(n_steps):
            # Update Hx, Hy
            self.hx[:, :-1] = chx_a[:, :-1] * self.hx[:, :-1] - chx_b[:, :-1] * (
                self.ez[:, 1:] - self.ez[:, :-1]
            )
            self.hy[:-1, :] = chy_a[:-1, :] * self.hy[:-1, :] + chy_b[:-1, :] * (
                self.ez[1:, :] - self.ez[:-1, :]
            )

            # Update Split Ez
            self.ezx[1:, :] = cx_a[1:, :] * self.ezx[1:, :] + cx_b[1:, :] * (
                self.hy[1:, :] - self.hy[:-1, :]
            )
            self.ezy[:, 1:] = cy_a[:, 1:] * self.ezy[:, 1:] - cy_b[:, 1:] * (
                self.hx[:, 1:] - self.hx[:, :-1]
            )
            self.ez = self.ezx + self.ezy

    def field_energy(self) -> float:
        """Total EM energy in the grid."""
        return float(np.sum(self.ez**2) + np.sum(self.hx**2) + np.sum(self.hy**2))

    def field_at_point(self, x: int, y: int) -> float:
        """Read Ez at a grid point."""
        return float(self.ez[x, y])

    def cross_section(self, x: int) -> np.ndarray:
        """Return Ez cross-section at a given x position."""
        return self.ez[x, :].copy()

    def snapshot(self) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
        """Return copies of all field components."""
        return self.ez.copy(), self.hx.copy(), self.hy.copy()

set_waveguide(y_center, width_cells, refractive_index=3.48, x_start=0, x_end=None)

Define a horizontal waveguide stripe.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
def set_waveguide(
    self,
    y_center: int,
    width_cells: int,
    refractive_index: float = 3.48,
    x_start: int = 0,
    x_end: Optional[int] = None,
) -> None:
    """Define a horizontal waveguide stripe."""
    if refractive_index < 1.0:
        raise ValueError(f"Invalid refractive index: {refractive_index}. Must be >= 1.0.")

    x_end = x_end or self.nx
    x_start = max(0, min(self.nx, x_start))
    x_end = max(0, min(self.nx, x_end))
    y_lo = max(0, min(self.ny, y_center - width_cells // 2))
    y_hi = max(0, min(self.ny, y_center + width_cells // 2))

    self.n_map[x_start:x_end, y_lo:y_hi] = refractive_index

inject_source(x, y, wavelength_nm=1550.0, amplitude=1.0, sigma_cells=10)

Inject a 2D Gaussian source.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
def inject_source(
    self,
    x: int,
    y: int,
    wavelength_nm: float = 1550.0,
    amplitude: float = 1.0,
    sigma_cells: int = 10,
) -> None:
    """Inject a 2D Gaussian source."""
    if wavelength_nm <= 0.0:
        raise ValueError(f"Wavelength must be strictly positive, got {wavelength_nm}")
    if not (0 <= x < self.nx) or not (0 <= y < self.ny):
        raise ValueError(f"Source injection ({x}, {y}) out of bounds [{self.nx}, {self.ny}]")

    freq = self.c0 / (wavelength_nm * 1e-9)
    for ix in range(max(0, x - 3 * sigma_cells), min(self.nx, x + 3 * sigma_cells)):
        for iy in range(max(0, y - 3 * sigma_cells), min(self.ny, y + 3 * sigma_cells)):
            dx_r = (ix - x) / sigma_cells
            dy_r = (iy - y) / sigma_cells
            envelope = amplitude * math.exp(-0.5 * (dx_r**2 + dy_r**2))
            self.ez[ix, iy] = envelope * math.cos(2 * math.pi * freq * 0)

step(n_steps=1)

Advance TE simulation by n_steps.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
def step(self, n_steps: int = 1) -> None:
    """Advance TE simulation by n_steps."""
    if np.any(self.n_map <= 0):
        raise ValueError("Refractive index must be > 0 in all cells.")

    eps0 = 8.854e-12
    mu0 = 4 * math.pi * 1e-7

    eps_map = eps0 * self.n_map**2

    cx_a = (eps_map - self.sigma_x * self.dt / 2.0) / (eps_map + self.sigma_x * self.dt / 2.0)
    cx_b = self.dt / ((eps_map + self.sigma_x * self.dt / 2.0) * self.dx)
    cy_a = (eps_map - self.sigma_y * self.dt / 2.0) / (eps_map + self.sigma_y * self.dt / 2.0)
    cy_b = self.dt / ((eps_map + self.sigma_y * self.dt / 2.0) * self.dy)

    smag_y = self.sigma_y * (mu0 / eps0)
    smag_x = self.sigma_x * (mu0 / eps0)
    chx_a = (mu0 - smag_y * self.dt / 2.0) / (mu0 + smag_y * self.dt / 2.0)
    chx_b = self.dt / ((mu0 + smag_y * self.dt / 2.0) * self.dy)
    chy_a = (mu0 - smag_x * self.dt / 2.0) / (mu0 + smag_x * self.dt / 2.0)
    chy_b = self.dt / ((mu0 + smag_x * self.dt / 2.0) * self.dx)

    for _ in range(n_steps):
        # Update Hx, Hy
        self.hx[:, :-1] = chx_a[:, :-1] * self.hx[:, :-1] - chx_b[:, :-1] * (
            self.ez[:, 1:] - self.ez[:, :-1]
        )
        self.hy[:-1, :] = chy_a[:-1, :] * self.hy[:-1, :] + chy_b[:-1, :] * (
            self.ez[1:, :] - self.ez[:-1, :]
        )

        # Update Split Ez
        self.ezx[1:, :] = cx_a[1:, :] * self.ezx[1:, :] + cx_b[1:, :] * (
            self.hy[1:, :] - self.hy[:-1, :]
        )
        self.ezy[:, 1:] = cy_a[:, 1:] * self.ezy[:, 1:] - cy_b[:, 1:] * (
            self.hx[:, 1:] - self.hx[:, :-1]
        )
        self.ez = self.ezx + self.ezy

field_energy()

Total EM energy in the grid.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
717
718
719
def field_energy(self) -> float:
    """Total EM energy in the grid."""
    return float(np.sum(self.ez**2) + np.sum(self.hx**2) + np.sum(self.hy**2))

field_at_point(x, y)

Read Ez at a grid point.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
721
722
723
def field_at_point(self, x: int, y: int) -> float:
    """Read Ez at a grid point."""
    return float(self.ez[x, y])

cross_section(x)

Return Ez cross-section at a given x position.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
725
726
727
def cross_section(self, x: int) -> np.ndarray:
    """Return Ez cross-section at a given x position."""
    return self.ez[x, :].copy()

snapshot()

Return copies of all field components.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
729
730
731
def snapshot(self) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Return copies of all field components."""
    return self.ez.copy(), self.hx.copy(), self.hy.copy()

MeepAdapter

Optional adapter for the Meep FDTD solver.

Provides a thin wrapper that translates SC-NeuroCore photonic configurations into Meep simulation objects. Meep is an optional dependency; all methods gracefully degrade if not installed.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
class MeepAdapter:
    """Optional adapter for the Meep FDTD solver.

    Provides a thin wrapper that translates SC-NeuroCore photonic
    configurations into Meep simulation objects. Meep is an optional
    dependency; all methods gracefully degrade if not installed.
    """

    @staticmethod
    def is_available() -> bool:
        """Check if Meep is installed."""
        try:
            import meep  # noqa: F401

            return True
        except ImportError:
            return False

    @staticmethod
    def build_waveguide_geometry(
        target: PhotonicTarget,
        waveguide_width_um: float = 0.5,
        length_um: float = 10.0,
        substrate_index: float = 1.45,
    ) -> Dict[str, Any]:
        """Build a Meep geometry dict (usable even without Meep).

        Returns a serialisable description of the simulation setup.
        """
        core_index = 3.48 if target.wavelength_nm > 1000 else 2.0
        wavelength_um = target.wavelength_nm / 1000.0
        freq = 1.0 / wavelength_um  # Meep normalised frequency

        return {
            "cell_size": [length_um, 3.0 * waveguide_width_um, 0],
            "resolution": 20,
            "sources": [
                {
                    "type": "ContinuousSource"
                    if target.modulation == OpticalModulation.PHASE
                    else "GaussianSource",
                    "frequency": freq,
                    "center": [-length_um / 2 + 0.5, 0, 0],
                    "size": [0, waveguide_width_um, 0],
                }
            ],
            "geometry": [
                {
                    "type": "Block",
                    "material_index": core_index,
                    "center": [0, 0, 0],
                    "size": [length_um, waveguide_width_um, "Infinity"],
                },
            ],
            "substrate_index": substrate_index,
            "pml_layers": 1.0,
            "wavelength_nm": target.wavelength_nm,
            "modulation": target.modulation.value,
        }

    @staticmethod
    def run_simulation(geometry: Dict[str, Any], run_time: float = 50.0) -> Dict[str, Any]:
        """Run a Meep simulation or return a mock result if unavailable."""
        if not MeepAdapter.is_available():
            # Mock result for testing without Meep
            return {
                "transmission": 0.85,
                "reflection": 0.02,
                "field_decay": 1e-4,
                "run_time": run_time,
                "mock": True,
                "wavelength_nm": geometry.get("wavelength_nm", 1550.0),
            }

        import meep as mp

        cell_size = geometry["cell_size"]
        resolution = geometry["resolution"]
        src_spec = geometry["sources"][0]
        geo_spec = geometry["geometry"][0]

        cell = mp.Vector3(*cell_size)
        freq = src_spec["frequency"]

        sources = [
            mp.Source(
                mp.ContinuousSource(frequency=freq),
                component=mp.Ez,
                center=mp.Vector3(*src_spec["center"]),
                size=mp.Vector3(*src_spec["size"]),
            )
        ]

        mat = mp.Medium(index=geo_spec["material_index"])
        geo_objs = [
            mp.Block(
                size=mp.Vector3(geo_spec["size"][0], geo_spec["size"][1]),
                center=mp.Vector3(*geo_spec["center"]),
                material=mat,
            )
        ]

        sim = mp.Simulation(
            cell_size=cell,
            resolution=resolution,
            sources=sources,
            geometry=geo_objs,
            boundary_layers=[mp.PML(geometry["pml_layers"])],
        )

        flux_region = mp.FluxRegion(
            center=mp.Vector3(cell_size[0] / 2 - 1, 0),
            size=mp.Vector3(0, cell_size[1]),
        )
        trans = sim.add_flux(freq, 0, 1, flux_region)

        sim.run(until=run_time)

        flux_data = mp.get_fluxes(trans)
        return {
            "transmission": float(flux_data[0]) if flux_data else 0.0,
            "reflection": 0.0,
            "field_decay": 0.0,
            "run_time": run_time,
            "mock": False,
            "wavelength_nm": geometry.get("wavelength_nm", 1550.0),
        }

is_available() staticmethod

Check if Meep is installed.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
745
746
747
748
749
750
751
752
753
@staticmethod
def is_available() -> bool:
    """Check if Meep is installed."""
    try:
        import meep  # noqa: F401

        return True
    except ImportError:
        return False

build_waveguide_geometry(target, waveguide_width_um=0.5, length_um=10.0, substrate_index=1.45) staticmethod

Build a Meep geometry dict (usable even without Meep).

Returns a serialisable description of the simulation setup.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
@staticmethod
def build_waveguide_geometry(
    target: PhotonicTarget,
    waveguide_width_um: float = 0.5,
    length_um: float = 10.0,
    substrate_index: float = 1.45,
) -> Dict[str, Any]:
    """Build a Meep geometry dict (usable even without Meep).

    Returns a serialisable description of the simulation setup.
    """
    core_index = 3.48 if target.wavelength_nm > 1000 else 2.0
    wavelength_um = target.wavelength_nm / 1000.0
    freq = 1.0 / wavelength_um  # Meep normalised frequency

    return {
        "cell_size": [length_um, 3.0 * waveguide_width_um, 0],
        "resolution": 20,
        "sources": [
            {
                "type": "ContinuousSource"
                if target.modulation == OpticalModulation.PHASE
                else "GaussianSource",
                "frequency": freq,
                "center": [-length_um / 2 + 0.5, 0, 0],
                "size": [0, waveguide_width_um, 0],
            }
        ],
        "geometry": [
            {
                "type": "Block",
                "material_index": core_index,
                "center": [0, 0, 0],
                "size": [length_um, waveguide_width_um, "Infinity"],
            },
        ],
        "substrate_index": substrate_index,
        "pml_layers": 1.0,
        "wavelength_nm": target.wavelength_nm,
        "modulation": target.modulation.value,
    }

run_simulation(geometry, run_time=50.0) staticmethod

Run a Meep simulation or return a mock result if unavailable.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
@staticmethod
def run_simulation(geometry: Dict[str, Any], run_time: float = 50.0) -> Dict[str, Any]:
    """Run a Meep simulation or return a mock result if unavailable."""
    if not MeepAdapter.is_available():
        # Mock result for testing without Meep
        return {
            "transmission": 0.85,
            "reflection": 0.02,
            "field_decay": 1e-4,
            "run_time": run_time,
            "mock": True,
            "wavelength_nm": geometry.get("wavelength_nm", 1550.0),
        }

    import meep as mp

    cell_size = geometry["cell_size"]
    resolution = geometry["resolution"]
    src_spec = geometry["sources"][0]
    geo_spec = geometry["geometry"][0]

    cell = mp.Vector3(*cell_size)
    freq = src_spec["frequency"]

    sources = [
        mp.Source(
            mp.ContinuousSource(frequency=freq),
            component=mp.Ez,
            center=mp.Vector3(*src_spec["center"]),
            size=mp.Vector3(*src_spec["size"]),
        )
    ]

    mat = mp.Medium(index=geo_spec["material_index"])
    geo_objs = [
        mp.Block(
            size=mp.Vector3(geo_spec["size"][0], geo_spec["size"][1]),
            center=mp.Vector3(*geo_spec["center"]),
            material=mat,
        )
    ]

    sim = mp.Simulation(
        cell_size=cell,
        resolution=resolution,
        sources=sources,
        geometry=geo_objs,
        boundary_layers=[mp.PML(geometry["pml_layers"])],
    )

    flux_region = mp.FluxRegion(
        center=mp.Vector3(cell_size[0] / 2 - 1, 0),
        size=mp.Vector3(0, cell_size[1]),
    )
    trans = sim.add_flux(freq, 0, 1, flux_region)

    sim.run(until=run_time)

    flux_data = mp.get_fluxes(trans)
    return {
        "transmission": float(flux_data[0]) if flux_data else 0.0,
        "reflection": 0.0,
        "field_decay": 0.0,
        "run_time": run_time,
        "mock": False,
        "wavelength_nm": geometry.get("wavelength_nm", 1550.0),
    }

WaveguidePair dataclass

A pair of adjacent waveguides for crosstalk analysis.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
@dataclass
class WaveguidePair:
    """A pair of adjacent waveguides for crosstalk analysis."""

    waveguide_width_nm: float = 450.0
    gap_nm: float = 200.0
    coupling_length_um: float = 10.0
    core_index: float = 3.48
    cladding_index: float = 1.45
    wavelength_nm: float = 1550.0

    @property
    def effective_index_diff(self) -> float:
        """Approximate effective index difference between even/odd modes."""
        # Exponential evanescent decay model
        decay_length_nm = self.wavelength_nm / (
            2 * math.pi * math.sqrt(self.core_index**2 - self.cladding_index**2)
        )
        return 0.1 * math.exp(-self.gap_nm / decay_length_nm)

    @property
    def coupling_coefficient(self) -> float:
        """Coupling coefficient κ (per micrometre)."""
        dn = self.effective_index_diff
        return math.pi * dn / (self.wavelength_nm * 1e-3)

    @property
    def coupling_ratio(self) -> float:
        """Power coupling ratio at the end of the coupling length."""
        kl = self.coupling_coefficient * self.coupling_length_um
        return math.sin(kl) ** 2

    @property
    def isolation_db(self) -> float:
        """Isolation in dB (lower coupling = better isolation)."""
        ratio = self.coupling_ratio
        if ratio < 1e-15:
            return 300.0
        return -10.0 * math.log10(max(ratio, 1e-30))

effective_index_diff property

Approximate effective index difference between even/odd modes.

coupling_coefficient property

Coupling coefficient κ (per micrometre).

coupling_ratio property

Power coupling ratio at the end of the coupling length.

isolation_db property

Isolation in dB (lower coupling = better isolation).

CrosstalkModel

Models evanescent crosstalk between adjacent waveguides.

Uses coupled-mode theory to compute power transfer between parallel waveguide runs on a photonic chip.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
 910
 911
 912
 913
 914
 915
 916
 917
 918
 919
 920
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
class CrosstalkModel:
    """Models evanescent crosstalk between adjacent waveguides.

    Uses coupled-mode theory to compute power transfer between
    parallel waveguide runs on a photonic chip.
    """

    def __init__(self) -> None:
        self.pairs: List[WaveguidePair] = []

    def add_pair(self, pair: WaveguidePair) -> None:
        self.pairs.append(pair)

    def transfer_matrix(self, pair: WaveguidePair) -> np.ndarray:
        """2×2 transfer matrix for a directional coupler."""
        kl = pair.coupling_coefficient * pair.coupling_length_um
        c = math.cos(kl)
        s = math.sin(kl)
        return np.array([[c, 1j * s], [1j * s, c]])

    def compute_crosstalk(
        self, pair: WaveguidePair, input_power: Tuple[float, float] = (1.0, 0.0)
    ) -> Tuple[float, float]:
        """Compute output power on both waveguides."""
        t = self.transfer_matrix(pair)
        inp = np.array(input_power, dtype=complex)
        out = t @ inp
        return float(np.abs(out[0]) ** 2), float(np.abs(out[1]) ** 2)

    def worst_case_isolation(self) -> float:
        """Worst-case isolation across all pairs (minimum dB)."""
        if not self.pairs:
            return float("inf")
        return min(p.isolation_db for p in self.pairs)

    def analyze_bank(
        self,
        waveguides: int,
        gap_nm: float,
        coupling_length_um: float,
        wavelength_nm: float = 1550.0,
        core_index: float = 3.48,
        cladding_index: float = 1.45,
    ) -> Dict[str, Any]:
        """Geometric crosstalk for a uniform bank of parallel waveguides.

        Uses coupled-mode theory with a Marcatili-form evanescent decay:

        ``L_decay = λ / (2π √(n_core² − n_clad²))``, ``κ = π Δn_eff / λ``,
        ``ratio = sin²(κL)``, ``isolation [dB] = −10 log₁₀(ratio)``.

        Adjacent pairs (gap = ``gap_nm``) are the dominant term; next-nearest
        pairs (gap = ``2·gap_nm``) are included as the largest secondary
        term. All further pairs decay at least as ``exp(−2·g/L_decay)``.

        Delegated to the Rust engine (``py_ph_analyze_crosstalk_bank``) when
        available; a Python fallback using :class:`WaveguidePair` matches
        the Rust result to within 1e-9.

        References:

        - Marcatili, Bell Syst. Tech. J. 48(7):2071-2102, 1969.
        - Okamoto, *Fundamentals of Optical Waveguides*, 2006, Ch. 4.
        """
        if waveguides < 1:
            raise ValueError("waveguides must be >= 1")

        if _HAS_RUST_PH:
            return dict(
                py_ph_analyze_crosstalk_bank(
                    num_waveguides=waveguides,
                    gap_nm=gap_nm,
                    coupling_length_um=coupling_length_um,
                    wavelength_nm=wavelength_nm,
                    core_index=core_index,
                    cladding_index=cladding_index,
                )
            )

        # Pure-Python fallback: identical math via WaveguidePair properties.
        near = WaveguidePair(
            gap_nm=gap_nm,
            coupling_length_um=coupling_length_um,
            wavelength_nm=wavelength_nm,
            core_index=core_index,
            cladding_index=cladding_index,
        )
        far = WaveguidePair(
            gap_nm=2.0 * gap_nm,
            coupling_length_um=coupling_length_um,
            wavelength_nm=wavelength_nm,
            core_index=core_index,
            cladding_index=cladding_index,
        )
        num_near = max(0, waveguides - 1)
        num_far = max(0, waveguides - 2)
        total = num_near + num_far
        if total == 0:
            worst = float("inf")
            mean_ratio = 0.0
            max_ratio = 0.0
        else:
            worst = min(near.isolation_db, far.isolation_db)
            mean_ratio = (num_near * near.coupling_ratio + num_far * far.coupling_ratio) / total
            max_ratio = max(near.coupling_ratio, far.coupling_ratio)
        return {
            "num_waveguides": waveguides,
            "num_pairs": num_near + num_far,
            "num_near_pairs": num_near,
            "num_far_pairs": num_far,
            "gap_nm": gap_nm,
            "coupling_length_um": coupling_length_um,
            "adjacent_coupling_ratio": near.coupling_ratio,
            "adjacent_isolation_db": near.isolation_db,
            "next_nearest_coupling_ratio": far.coupling_ratio,
            "next_nearest_isolation_db": far.isolation_db,
            "worst_isolation_db": worst,
            "mean_coupling_ratio": mean_ratio,
            "max_coupling_ratio": max_ratio,
            "crosstalk_safe": worst > 20.0,
            "backend": "python",
        }

    def analyze_pairs(
        self,
        pair_indices: List[Tuple[int, int]],
        gaps_nm: List[float],
        coupling_lengths_um: List[float],
        wavelength_nm: float = 1550.0,
        core_index: float = 3.48,
        cladding_index: float = 1.45,
    ) -> Dict[str, Any]:
        """Per-pair crosstalk for arbitrary waveguide geometry — the O(N²) path.

        Use this when the bank is not uniform: each pair may have its own
        gap and coupling length. Rust path evaluates pairs in parallel via
        Rayon; the Python fallback reproduces the same math serially.
        """
        n = len(pair_indices)
        if len(gaps_nm) != n or len(coupling_lengths_um) != n:
            raise ValueError(
                f"pair_indices ({n}), gaps_nm ({len(gaps_nm)}) and "
                f"coupling_lengths_um ({len(coupling_lengths_um)}) must be equal length"
            )

        if _HAS_RUST_PH and n > 0:
            pairs_a = [a for (a, _) in pair_indices]
            pairs_b = [b for (_, b) in pair_indices]
            return dict(
                py_ph_analyze_crosstalk_pairs(
                    pairs_a=pairs_a,
                    pairs_b=pairs_b,
                    gaps_nm=list(gaps_nm),
                    lengths_um=list(coupling_lengths_um),
                    wavelength_nm=wavelength_nm,
                    core_index=core_index,
                    cladding_index=cladding_index,
                )
            )

        # Pure-Python fallback.
        kappas: List[float] = []
        ratios: List[float] = []
        isos: List[float] = []
        for g, ell in zip(gaps_nm, coupling_lengths_um):
            p = WaveguidePair(
                gap_nm=g,
                coupling_length_um=ell,
                wavelength_nm=wavelength_nm,
                core_index=core_index,
                cladding_index=cladding_index,
            )
            kappas.append(p.coupling_coefficient)
            ratios.append(p.coupling_ratio)
            isos.append(p.isolation_db)
        return {
            "pair_a": [a for (a, _) in pair_indices],
            "pair_b": [b for (_, b) in pair_indices],
            "gap_nm": list(gaps_nm),
            "coupling_length_um": list(coupling_lengths_um),
            "coupling_coefficient_per_um": kappas,
            "coupling_ratio": ratios,
            "isolation_db": isos,
            "num_pairs": n,
            "backend": "python",
        }

transfer_matrix(pair)

2×2 transfer matrix for a directional coupler.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
923
924
925
926
927
928
def transfer_matrix(self, pair: WaveguidePair) -> np.ndarray:
    """2×2 transfer matrix for a directional coupler."""
    kl = pair.coupling_coefficient * pair.coupling_length_um
    c = math.cos(kl)
    s = math.sin(kl)
    return np.array([[c, 1j * s], [1j * s, c]])

compute_crosstalk(pair, input_power=(1.0, 0.0))

Compute output power on both waveguides.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
930
931
932
933
934
935
936
937
def compute_crosstalk(
    self, pair: WaveguidePair, input_power: Tuple[float, float] = (1.0, 0.0)
) -> Tuple[float, float]:
    """Compute output power on both waveguides."""
    t = self.transfer_matrix(pair)
    inp = np.array(input_power, dtype=complex)
    out = t @ inp
    return float(np.abs(out[0]) ** 2), float(np.abs(out[1]) ** 2)

worst_case_isolation()

Worst-case isolation across all pairs (minimum dB).

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
939
940
941
942
943
def worst_case_isolation(self) -> float:
    """Worst-case isolation across all pairs (minimum dB)."""
    if not self.pairs:
        return float("inf")
    return min(p.isolation_db for p in self.pairs)

analyze_bank(waveguides, gap_nm, coupling_length_um, wavelength_nm=1550.0, core_index=3.48, cladding_index=1.45)

Geometric crosstalk for a uniform bank of parallel waveguides.

Uses coupled-mode theory with a Marcatili-form evanescent decay:

L_decay = λ / (2π √(n_core² − n_clad²)), κ = π Δn_eff / λ, ratio = sin²(κL), isolation [dB] = −10 log₁₀(ratio).

Adjacent pairs (gap = gap_nm) are the dominant term; next-nearest pairs (gap = 2·gap_nm) are included as the largest secondary term. All further pairs decay at least as exp(−2·g/L_decay).

Delegated to the Rust engine (py_ph_analyze_crosstalk_bank) when available; a Python fallback using :class:WaveguidePair matches the Rust result to within 1e-9.

References:

  • Marcatili, Bell Syst. Tech. J. 48(7):2071-2102, 1969.
  • Okamoto, Fundamentals of Optical Waveguides, 2006, Ch. 4.
Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
def analyze_bank(
    self,
    waveguides: int,
    gap_nm: float,
    coupling_length_um: float,
    wavelength_nm: float = 1550.0,
    core_index: float = 3.48,
    cladding_index: float = 1.45,
) -> Dict[str, Any]:
    """Geometric crosstalk for a uniform bank of parallel waveguides.

    Uses coupled-mode theory with a Marcatili-form evanescent decay:

    ``L_decay = λ / (2π √(n_core² − n_clad²))``, ``κ = π Δn_eff / λ``,
    ``ratio = sin²(κL)``, ``isolation [dB] = −10 log₁₀(ratio)``.

    Adjacent pairs (gap = ``gap_nm``) are the dominant term; next-nearest
    pairs (gap = ``2·gap_nm``) are included as the largest secondary
    term. All further pairs decay at least as ``exp(−2·g/L_decay)``.

    Delegated to the Rust engine (``py_ph_analyze_crosstalk_bank``) when
    available; a Python fallback using :class:`WaveguidePair` matches
    the Rust result to within 1e-9.

    References:

    - Marcatili, Bell Syst. Tech. J. 48(7):2071-2102, 1969.
    - Okamoto, *Fundamentals of Optical Waveguides*, 2006, Ch. 4.
    """
    if waveguides < 1:
        raise ValueError("waveguides must be >= 1")

    if _HAS_RUST_PH:
        return dict(
            py_ph_analyze_crosstalk_bank(
                num_waveguides=waveguides,
                gap_nm=gap_nm,
                coupling_length_um=coupling_length_um,
                wavelength_nm=wavelength_nm,
                core_index=core_index,
                cladding_index=cladding_index,
            )
        )

    # Pure-Python fallback: identical math via WaveguidePair properties.
    near = WaveguidePair(
        gap_nm=gap_nm,
        coupling_length_um=coupling_length_um,
        wavelength_nm=wavelength_nm,
        core_index=core_index,
        cladding_index=cladding_index,
    )
    far = WaveguidePair(
        gap_nm=2.0 * gap_nm,
        coupling_length_um=coupling_length_um,
        wavelength_nm=wavelength_nm,
        core_index=core_index,
        cladding_index=cladding_index,
    )
    num_near = max(0, waveguides - 1)
    num_far = max(0, waveguides - 2)
    total = num_near + num_far
    if total == 0:
        worst = float("inf")
        mean_ratio = 0.0
        max_ratio = 0.0
    else:
        worst = min(near.isolation_db, far.isolation_db)
        mean_ratio = (num_near * near.coupling_ratio + num_far * far.coupling_ratio) / total
        max_ratio = max(near.coupling_ratio, far.coupling_ratio)
    return {
        "num_waveguides": waveguides,
        "num_pairs": num_near + num_far,
        "num_near_pairs": num_near,
        "num_far_pairs": num_far,
        "gap_nm": gap_nm,
        "coupling_length_um": coupling_length_um,
        "adjacent_coupling_ratio": near.coupling_ratio,
        "adjacent_isolation_db": near.isolation_db,
        "next_nearest_coupling_ratio": far.coupling_ratio,
        "next_nearest_isolation_db": far.isolation_db,
        "worst_isolation_db": worst,
        "mean_coupling_ratio": mean_ratio,
        "max_coupling_ratio": max_ratio,
        "crosstalk_safe": worst > 20.0,
        "backend": "python",
    }

analyze_pairs(pair_indices, gaps_nm, coupling_lengths_um, wavelength_nm=1550.0, core_index=3.48, cladding_index=1.45)

Per-pair crosstalk for arbitrary waveguide geometry — the O(N²) path.

Use this when the bank is not uniform: each pair may have its own gap and coupling length. Rust path evaluates pairs in parallel via Rayon; the Python fallback reproduces the same math serially.

Source code in src/sc_neurocore/optics/photonic_emitter.py
Python
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
def analyze_pairs(
    self,
    pair_indices: List[Tuple[int, int]],
    gaps_nm: List[float],
    coupling_lengths_um: List[float],
    wavelength_nm: float = 1550.0,
    core_index: float = 3.48,
    cladding_index: float = 1.45,
) -> Dict[str, Any]:
    """Per-pair crosstalk for arbitrary waveguide geometry — the O(N²) path.

    Use this when the bank is not uniform: each pair may have its own
    gap and coupling length. Rust path evaluates pairs in parallel via
    Rayon; the Python fallback reproduces the same math serially.
    """
    n = len(pair_indices)
    if len(gaps_nm) != n or len(coupling_lengths_um) != n:
        raise ValueError(
            f"pair_indices ({n}), gaps_nm ({len(gaps_nm)}) and "
            f"coupling_lengths_um ({len(coupling_lengths_um)}) must be equal length"
        )

    if _HAS_RUST_PH and n > 0:
        pairs_a = [a for (a, _) in pair_indices]
        pairs_b = [b for (_, b) in pair_indices]
        return dict(
            py_ph_analyze_crosstalk_pairs(
                pairs_a=pairs_a,
                pairs_b=pairs_b,
                gaps_nm=list(gaps_nm),
                lengths_um=list(coupling_lengths_um),
                wavelength_nm=wavelength_nm,
                core_index=core_index,
                cladding_index=cladding_index,
            )
        )

    # Pure-Python fallback.
    kappas: List[float] = []
    ratios: List[float] = []
    isos: List[float] = []
    for g, ell in zip(gaps_nm, coupling_lengths_um):
        p = WaveguidePair(
            gap_nm=g,
            coupling_length_um=ell,
            wavelength_nm=wavelength_nm,
            core_index=core_index,
            cladding_index=cladding_index,
        )
        kappas.append(p.coupling_coefficient)
        ratios.append(p.coupling_ratio)
        isos.append(p.isolation_db)
    return {
        "pair_a": [a for (a, _) in pair_indices],
        "pair_b": [b for (_, b) in pair_indices],
        "gap_nm": list(gaps_nm),
        "coupling_length_um": list(coupling_lengths_um),
        "coupling_coefficient_per_um": kappas,
        "coupling_ratio": ratios,
        "isolation_db": isos,
        "num_pairs": n,
        "backend": "python",
    }