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. The stable photonic_emitter import path is a definition-free compatibility facade; seven bounded modules own types, conversion, FDTD, compilation, Meep, crosstalk, and emitter responsibilities.

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 crosstalk analysis; a validated pure-Python fallback uses the same closed-form model when the engine wheel is absent. Standalone Rust, Go, Julia, and Mojo implementations mirror only that numeric crosstalk contract. FDTD, Meep, compilation, netlist emission, and filesystem work remain owned by Python and are not presented as cross-language capabilities.

---

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; unavailable runtimes raise ImportError and never fabricate simulation results
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
Native crosstalk parity	Engine Rust and standalone Rust/Go/Julia/Mojo execute source-bound parity contracts against Python
Focused verification	102 tests; exact 100% coverage over 702 statements and 214 branches in the facade and seven responsibility modules

---

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

The validated Python fallback is exercised when _HAS_RUST_PH == False. Focused engine tests require relative tolerance 1e-9 and absolute tolerance 1e-12. The separate source-bound Rust, Go, Julia, and Mojo crosstalk mirrors compile and execute against tighter per-runtime envelopes; they do not mirror the routing, MZI, power-budget, FDTD, Meep, netlist, or GDS surfaces.

---

7. Performance benchmarks

The committed crosstalk harness is benchmarks/bench_photonic_crosstalk.py. It validates source hashes and the first output tuple before recording timing data in benchmarks/results/local_python_2026-07-14_photonic_crosstalk.json.

7.1 Enrolled workload and parity

The workload contains 4,096 pairs. Gap is 180 + index mod 64 nm and coupling length is 8 + index mod 17 µm at 1550 nm with core/cladding indices 3.48/1.45. The final measured first-pair maximum absolute differences from Python were:

Runtime	Maximum absolute difference	Enrolled envelope
Rust	0	1e-15
Go	1.78e-15	3e-15
Julia	0	1e-15
Mojo	1.35e-10	2e-10

The runtime-specific envelopes reflect ordinary floating-point-library and ABI differences. They are correctness bounds, not a claim that the implementations are bit-identical.

7.2 Local regression timings

Runtime	Median time per 4,096-pair batch
Python	12.785 ms
Rust	0.165 ms
Go	0.422 ms
Julia	0.219 ms
Mojo	0.144 ms

These results were recorded on Linux x86-64 with CPython 3.12.3, Rust 1.96.0, Go 1.24.0, Julia 1.12.6, and Mojo 1.0.0b1. The process was pinned to logical CPU 0, but the host had no isolated CPUs, used the powersave governor, and was loaded. Native runtimes also use different call boundaries. The artefact is therefore marked local_regression_non_isolated and promotion_eligible: false: it is reproducible local diagnostic evidence, not a universal speed or production-throughput claim.

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]. The optional-extra CI job installs and executes the GDSII tests; minimal environments skip that dependency-gated slice.
• MeepAdapter.run_simulation requires pymeep. An unavailable Meep runtime raises ImportError. The adapter does not return synthetic field, flux, or energy data as a substitute for a simulation.
• 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

• Stable facade: src/sc_neurocore/optics/photonic_emitter.py (64 lines).
• Responsibility modules: src/sc_neurocore/optics/_photonic_{types,conversion,fdtd,compiler,meep,crosstalk,emitter}.py (106, 117, 297, 281, 218, 321, and 83 lines respectively). Their dependency graph and per-module ceilings are executable architecture contracts.
• Layer module: src/sc_neurocore/optics/photonic_layer.py.
• Focused verification: 102 tests, including public compatibility, deterministic valid-output digests, validation, FDTD, GDSII, bridge integration, native engine parity, and executable standalone Rust/Go/Julia/Mojo crosstalk parity. The facade and seven responsibility modules have exact 100% statement and branch coverage over 702 statements and 214 branches.
• Rust engine: engine/src/photonic.rs; the crosstalk floor and zero-length contract are shared with the standalone mirrors.
• Benchmark: benchmarks/bench_photonic_crosstalk.py, with native-build helpers in benchmarks/_photonic_crosstalk_native.py and committed local evidence in benchmarks/results/local_python_2026-07-14_photonic_crosstalk.json.
• 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

@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

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

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

Stable photonic compiler, co-simulation, layout, and crosstalk API.

Implementation ownership is split by responsibility across private sibling modules. Historical imports, class identities, and pickle-qualified paths remain anchored here.

CompilationResult dataclass

Result of a photonic compilation pass.

Source code in src/sc_neurocore/optics/_photonic_compiler.py

@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 __post_init__(self) -> None:
        """Validate compilation metadata before export."""
        if not isinstance(self.target, str) or not self.target.strip():
            raise ValueError("target must be a non-empty string")
        _require_count(self.num_modulators, "num_modulators")
        _require_non_negative(self.optical_power_mean_mw, "optical_power_mean_mw")
        _require_non_negative(self.phase_coverage_rad, "phase_coverage_rad")
        _require_non_negative(self.fdtd_energy, "fdtd_energy")
        if not isinstance(self.netlist, str):
            raise TypeError("netlist must be a string")

    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 GDSII through gdsfactory.

        The layout is a linear cascade at ``pitch_um`` spacing. An identifying
        header and a bounded copy of the logical netlist are stored on TEXT
        layer 63/0. The optional ``gdsfactory`` dependency is loaded only when
        this method is called.
        """
        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)."
            )
        if not isinstance(filename, str) or not filename.strip():
            raise ValueError("filename must be a non-empty string")
        _require_positive(mzi_length_um, "mzi_length_um")
        _require_positive(pitch_um, "pitch_um")
        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

        try:
            gf.get_active_pdk()
        except (ValueError, AttributeError):
            try:
                gf.gpdk.PDK.activate()
            except AttributeError:
                from gdsfactory.gpdk import get_generic_pdk

                get_generic_pdk().activate()

        kdb_cell = gf.kcl.create_cell(f"SC_NeuroCore_Target_{self.target}", allow_duplicate=True)
        component = gf.Component(kdb_cell=kdb_cell)
        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):
            reference = component.add_ref(mzi_cell)
            reference.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,
        }

__post_init__()

Validate compilation metadata before export.

Source code in src/sc_neurocore/optics/_photonic_compiler.py

def __post_init__(self) -> None:
    """Validate compilation metadata before export."""
    if not isinstance(self.target, str) or not self.target.strip():
        raise ValueError("target must be a non-empty string")
    _require_count(self.num_modulators, "num_modulators")
    _require_non_negative(self.optical_power_mean_mw, "optical_power_mean_mw")
    _require_non_negative(self.phase_coverage_rad, "phase_coverage_rad")
    _require_non_negative(self.fdtd_energy, "fdtd_energy")
    if not isinstance(self.netlist, str):
        raise TypeError("netlist must be a string")

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

Export the compiled MZI cascade to GDSII through gdsfactory.

The layout is a linear cascade at pitch_um spacing. An identifying header and a bounded copy of the logical netlist are stored on TEXT layer 63/0. The optional gdsfactory dependency is loaded only when this method is called.

Source code in src/sc_neurocore/optics/_photonic_compiler.py

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 GDSII through gdsfactory.

    The layout is a linear cascade at ``pitch_um`` spacing. An identifying
    header and a bounded copy of the logical netlist are stored on TEXT
    layer 63/0. The optional ``gdsfactory`` dependency is loaded only when
    this method is called.
    """
    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)."
        )
    if not isinstance(filename, str) or not filename.strip():
        raise ValueError("filename must be a non-empty string")
    _require_positive(mzi_length_um, "mzi_length_um")
    _require_positive(pitch_um, "pitch_um")
    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

    try:
        gf.get_active_pdk()
    except (ValueError, AttributeError):
        try:
            gf.gpdk.PDK.activate()
        except AttributeError:
            from gdsfactory.gpdk import get_generic_pdk

            get_generic_pdk().activate()

    kdb_cell = gf.kcl.create_cell(f"SC_NeuroCore_Target_{self.target}", allow_duplicate=True)
    component = gf.Component(kdb_cell=kdb_cell)
    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):
        reference = component.add_ref(mzi_cell)
        reference.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

Compile an SC bitstream into optical mapping, netlist, and co-simulation.

Source code in src/sc_neurocore/optics/_photonic_compiler.py

class PhotonicCompiler:
    """Compile an SC bitstream into optical mapping, netlist, and co-simulation."""

    def __init__(self, target: Optional[PhotonicTarget] = None):
        if target is not None and not isinstance(target, PhotonicTarget):
            raise TypeError("target must be a PhotonicTarget or 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[Any, Any],
        run_fdtd: bool = False,
        fdtd_steps: int = 100,
    ) -> CompilationResult:
        """Compile one non-empty binary SC bitstream to a photonic deployment."""
        if bitstream is None:
            raise ValueError("Input bitstream cannot be empty.")
        if not isinstance(run_fdtd, bool):
            raise TypeError("run_fdtd must be a Boolean")
        _require_count(fdtd_steps, "fdtd_steps")
        normalised = _normalise_bitstream(bitstream)
        if normalised.size == 0:
            raise ValueError("Input bitstream cannot be empty.")

        phases = self.converter.to_phase_array(normalised)
        power = self.converter.optical_power_profile(normalised)
        amplitudes = self.converter.to_amplitude_array(normalised)
        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 index, (phase, amplitude) in enumerate(zip(phases, amplitudes)):
            if self.target.modulator_type == "MZI":
                netlist_lines.append(f"ADD MZI mod_{index}")
                netlist_lines.append(f"SET mod_{index}:phase {phase:.6f}")
                netlist_lines.append(f"SET mod_{index}:amplitude {amplitude:.6f}")
            else:
                netlist_lines.append(f"ADD MICRORING ring_{index}")
                netlist_lines.append(f"SET ring_{index}:coupling {amplitude:.6f}")
                netlist_lines.append(f"SET ring_{index}:detuning {phase:.6f}")

        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="\n".join(netlist_lines),
            fdtd_energy=fdtd_energy,
        )

    def generate_mzi_verilog(self, bit_width: int = 16) -> str:
        """Generate SystemVerilog for an MZI modulator."""
        _require_count(bit_width, "bit_width", minimum=2)
        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."""
        _require_count(bit_width, "bit_width", minimum=2)
        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 one non-empty binary SC bitstream to a photonic deployment.

Source code in src/sc_neurocore/optics/_photonic_compiler.py

def compile_bitstream(
    self,
    bitstream: np.ndarray[Any, Any],
    run_fdtd: bool = False,
    fdtd_steps: int = 100,
) -> CompilationResult:
    """Compile one non-empty binary SC bitstream to a photonic deployment."""
    if bitstream is None:
        raise ValueError("Input bitstream cannot be empty.")
    if not isinstance(run_fdtd, bool):
        raise TypeError("run_fdtd must be a Boolean")
    _require_count(fdtd_steps, "fdtd_steps")
    normalised = _normalise_bitstream(bitstream)
    if normalised.size == 0:
        raise ValueError("Input bitstream cannot be empty.")

    phases = self.converter.to_phase_array(normalised)
    power = self.converter.optical_power_profile(normalised)
    amplitudes = self.converter.to_amplitude_array(normalised)
    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 index, (phase, amplitude) in enumerate(zip(phases, amplitudes)):
        if self.target.modulator_type == "MZI":
            netlist_lines.append(f"ADD MZI mod_{index}")
            netlist_lines.append(f"SET mod_{index}:phase {phase:.6f}")
            netlist_lines.append(f"SET mod_{index}:amplitude {amplitude:.6f}")
        else:
            netlist_lines.append(f"ADD MICRORING ring_{index}")
            netlist_lines.append(f"SET ring_{index}:coupling {amplitude:.6f}")
            netlist_lines.append(f"SET ring_{index}:detuning {phase:.6f}")

    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="\n".join(netlist_lines),
        fdtd_energy=fdtd_energy,
    )

generate_mzi_verilog(bit_width=16)

Generate SystemVerilog for an MZI modulator.

Source code in src/sc_neurocore/optics/_photonic_compiler.py

def generate_mzi_verilog(self, bit_width: int = 16) -> str:
    """Generate SystemVerilog for an MZI modulator."""
    _require_count(bit_width, "bit_width", minimum=2)
    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_compiler.py

def generate_microring_verilog(self, bit_width: int = 16) -> str:
    """Generate SystemVerilog for a microring resonator."""
    _require_count(bit_width, "bit_width", minimum=2)
    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
    """)

BitstreamToOptical

Convert SC bitstreams into optical pulse trains.

Source code in src/sc_neurocore/optics/_photonic_conversion.py

class BitstreamToOptical:
    """Convert SC bitstreams into optical pulse trains."""

    def __init__(self, target: PhotonicTarget):
        if not isinstance(target, PhotonicTarget):
            raise TypeError("target must be a PhotonicTarget")
        self.target = target

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

        Phase modulation maps one to phase zero and zero to phase π.
        Amplitude modulation maps one to unit amplitude and zero to zero.
        Hybrid modulation combines phase and amplitude encoding.
        """
        _require_positive(pulse_duration_ps, "pulse_duration_ps")
        pulses: List[OpticalPulse] = []
        for bit in _normalise_bitstream(bitstream):
            b = int(bit)
            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[Any, Any]) -> np.ndarray[Any, Any]:
        """Return the vectorised phase encoding in radians."""
        bs = _normalise_bitstream(bitstream)
        if self.target.modulation == OpticalModulation.PHASE:
            return np.where(bs > 0.5, 0.0, math.pi)
        if self.target.modulation == OpticalModulation.AMPLITUDE:
            return np.zeros_like(bs)
        return np.where(bs > 0.5, 0.0, math.pi / 2)

    def to_amplitude_array(self, bitstream: np.ndarray[Any, Any]) -> np.ndarray[Any, Any]:
        """Return the vectorised normalised amplitude encoding."""
        bs = _normalise_bitstream(bitstream)
        if self.target.modulation == OpticalModulation.PHASE:
            return np.ones_like(bs)
        if self.target.modulation == OpticalModulation.AMPLITUDE:
            return bs
        return 0.8 + 0.2 * bs

    def optical_power_profile(
        self,
        bitstream: np.ndarray[Any, Any],
        input_power_mw: float = 1.0,
    ) -> np.ndarray[Any, Any]:
        """Compute output power after the target insertion loss."""
        _require_non_negative(input_power_mw, "input_power_mw")
        amplitudes = self.to_amplitude_array(bitstream)
        loss_linear = 10.0 ** (-self.target.insertion_loss_db / 10.0)
        optical_power: np.ndarray[Any, Any] = amplitudes * amplitudes * input_power_mw * loss_linear
        return optical_power

convert(bitstream, pulse_duration_ps=10.0)

Map a binary SC bitstream to an optical pulse train.

Phase modulation maps one to phase zero and zero to phase π. Amplitude modulation maps one to unit amplitude and zero to zero. Hybrid modulation combines phase and amplitude encoding.

Source code in src/sc_neurocore/optics/_photonic_conversion.py

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

    Phase modulation maps one to phase zero and zero to phase π.
    Amplitude modulation maps one to unit amplitude and zero to zero.
    Hybrid modulation combines phase and amplitude encoding.
    """
    _require_positive(pulse_duration_ps, "pulse_duration_ps")
    pulses: List[OpticalPulse] = []
    for bit in _normalise_bitstream(bitstream):
        b = int(bit)
        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)

Return the vectorised phase encoding in radians.

Source code in src/sc_neurocore/optics/_photonic_conversion.py

def to_phase_array(self, bitstream: np.ndarray[Any, Any]) -> np.ndarray[Any, Any]:
    """Return the vectorised phase encoding in radians."""
    bs = _normalise_bitstream(bitstream)
    if self.target.modulation == OpticalModulation.PHASE:
        return np.where(bs > 0.5, 0.0, math.pi)
    if self.target.modulation == OpticalModulation.AMPLITUDE:
        return np.zeros_like(bs)
    return np.where(bs > 0.5, 0.0, math.pi / 2)

to_amplitude_array(bitstream)

Return the vectorised normalised amplitude encoding.

Source code in src/sc_neurocore/optics/_photonic_conversion.py

def to_amplitude_array(self, bitstream: np.ndarray[Any, Any]) -> np.ndarray[Any, Any]:
    """Return the vectorised normalised amplitude encoding."""
    bs = _normalise_bitstream(bitstream)
    if self.target.modulation == OpticalModulation.PHASE:
        return np.ones_like(bs)
    if self.target.modulation == OpticalModulation.AMPLITUDE:
        return bs
    return 0.8 + 0.2 * bs

optical_power_profile(bitstream, input_power_mw=1.0)

Compute output power after the target insertion loss.

Source code in src/sc_neurocore/optics/_photonic_conversion.py

def optical_power_profile(
    self,
    bitstream: np.ndarray[Any, Any],
    input_power_mw: float = 1.0,
) -> np.ndarray[Any, Any]:
    """Compute output power after the target insertion loss."""
    _require_non_negative(input_power_mw, "input_power_mw")
    amplitudes = self.to_amplitude_array(bitstream)
    loss_linear = 10.0 ** (-self.target.insertion_loss_db / 10.0)
    optical_power: np.ndarray[Any, Any] = amplitudes * amplitudes * input_power_mw * loss_linear
    return optical_power

CrosstalkModel

Evaluate evanescent crosstalk between parallel waveguide runs.

Source code in src/sc_neurocore/optics/_photonic_crosstalk.py

class CrosstalkModel:
    """Evaluate evanescent crosstalk between parallel waveguide runs."""

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

    def add_pair(self, pair: WaveguidePair) -> None:
        """Append one validated waveguide pair to the analyzer batch."""
        if not isinstance(pair, WaveguidePair):
            raise TypeError("pair must be a WaveguidePair")
        self.pairs.append(pair)

    def transfer_matrix(self, pair: WaveguidePair) -> np.ndarray[Any, Any]:
        """Return the two-by-two unitary directional-coupler matrix."""
        if not isinstance(pair, WaveguidePair):
            raise TypeError("pair must be a WaveguidePair")
        kl = pair.coupling_coefficient * pair.coupling_length_um
        cosine = math.cos(kl)
        sine = math.sin(kl)
        return np.array([[cosine, 1j * sine], [1j * sine, cosine]])

    def compute_crosstalk(
        self, pair: WaveguidePair, input_power: Tuple[float, float] = (1.0, 0.0)
    ) -> Tuple[float, float]:
        """Return output power on both waveguides for two input amplitudes."""
        if len(input_power) != 2:
            raise ValueError("input_power must contain exactly two field amplitudes")
        _require_finite(input_power[0], "input_power[0]")
        _require_finite(input_power[1], "input_power[1]")
        transfer = self.transfer_matrix(pair)
        output = transfer @ np.array(input_power, dtype=complex)
        return float(np.abs(output[0]) ** 2), float(np.abs(output[1]) ** 2)

    def worst_case_isolation(self) -> float:
        """Return minimum isolation across registered pairs in decibels."""
        if not self.pairs:
            return float("inf")
        return min(pair.isolation_db for pair 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]:
        """Analyze adjacent and next-nearest pairs in a uniform waveguide bank."""
        _require_index(waveguides, "waveguides")
        if waveguides < 1:
            raise ValueError("waveguides must be >= 1")
        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,
        )

        if _rust_backend_enabled():
            analyzer = _backend_analyzer(
                "py_ph_analyze_crosstalk_bank", py_ph_analyze_crosstalk_bank
            )
            return dict(
                analyzer(
                    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,
                )
            )

        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": total,
            "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]:
        """Analyze per-pair crosstalk for arbitrary waveguide geometry."""
        pair_count = len(pair_indices)
        if len(gaps_nm) != pair_count or len(coupling_lengths_um) != pair_count:
            raise ValueError(
                f"pair_indices ({pair_count}), gaps_nm ({len(gaps_nm)}) and "
                f"coupling_lengths_um ({len(coupling_lengths_um)}) must be equal length"
            )
        for index, (pair_a, pair_b) in enumerate(pair_indices):
            _require_index(pair_a, f"pair_indices[{index}][0]")
            _require_index(pair_b, f"pair_indices[{index}][1]")
            if pair_a == pair_b:
                raise ValueError(f"pair_indices[{index}] must name two distinct waveguides")

        pairs = [
            WaveguidePair(
                gap_nm=gap,
                coupling_length_um=length,
                wavelength_nm=wavelength_nm,
                core_index=core_index,
                cladding_index=cladding_index,
            )
            for gap, length in zip(gaps_nm, coupling_lengths_um)
        ]
        if _rust_backend_enabled() and pair_count > 0:
            analyzer = _backend_analyzer(
                "py_ph_analyze_crosstalk_pairs", py_ph_analyze_crosstalk_pairs
            )
            return dict(
                analyzer(
                    pairs_a=[pair_a for pair_a, _ in pair_indices],
                    pairs_b=[pair_b for _, pair_b in pair_indices],
                    gaps_nm=list(gaps_nm),
                    lengths_um=list(coupling_lengths_um),
                    wavelength_nm=wavelength_nm,
                    core_index=core_index,
                    cladding_index=cladding_index,
                )
            )

        return {
            "pair_a": [pair_a for pair_a, _ in pair_indices],
            "pair_b": [pair_b for _, pair_b in pair_indices],
            "gap_nm": list(gaps_nm),
            "coupling_length_um": list(coupling_lengths_um),
            "coupling_coefficient_per_um": [pair.coupling_coefficient for pair in pairs],
            "coupling_ratio": [pair.coupling_ratio for pair in pairs],
            "isolation_db": [pair.isolation_db for pair in pairs],
            "num_pairs": pair_count,
            "backend": "python",
        }

add_pair(pair)

Append one validated waveguide pair to the analyzer batch.

Source code in src/sc_neurocore/optics/_photonic_crosstalk.py

def add_pair(self, pair: WaveguidePair) -> None:
    """Append one validated waveguide pair to the analyzer batch."""
    if not isinstance(pair, WaveguidePair):
        raise TypeError("pair must be a WaveguidePair")
    self.pairs.append(pair)

transfer_matrix(pair)

Return the two-by-two unitary directional-coupler matrix.

Source code in src/sc_neurocore/optics/_photonic_crosstalk.py

def transfer_matrix(self, pair: WaveguidePair) -> np.ndarray[Any, Any]:
    """Return the two-by-two unitary directional-coupler matrix."""
    if not isinstance(pair, WaveguidePair):
        raise TypeError("pair must be a WaveguidePair")
    kl = pair.coupling_coefficient * pair.coupling_length_um
    cosine = math.cos(kl)
    sine = math.sin(kl)
    return np.array([[cosine, 1j * sine], [1j * sine, cosine]])

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

Return output power on both waveguides for two input amplitudes.

Source code in src/sc_neurocore/optics/_photonic_crosstalk.py

def compute_crosstalk(
    self, pair: WaveguidePair, input_power: Tuple[float, float] = (1.0, 0.0)
) -> Tuple[float, float]:
    """Return output power on both waveguides for two input amplitudes."""
    if len(input_power) != 2:
        raise ValueError("input_power must contain exactly two field amplitudes")
    _require_finite(input_power[0], "input_power[0]")
    _require_finite(input_power[1], "input_power[1]")
    transfer = self.transfer_matrix(pair)
    output = transfer @ np.array(input_power, dtype=complex)
    return float(np.abs(output[0]) ** 2), float(np.abs(output[1]) ** 2)

worst_case_isolation()

Return minimum isolation across registered pairs in decibels.

Source code in src/sc_neurocore/optics/_photonic_crosstalk.py

def worst_case_isolation(self) -> float:
    """Return minimum isolation across registered pairs in decibels."""
    if not self.pairs:
        return float("inf")
    return min(pair.isolation_db for pair in self.pairs)

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

Analyze adjacent and next-nearest pairs in a uniform waveguide bank.

Source code in src/sc_neurocore/optics/_photonic_crosstalk.py

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]:
    """Analyze adjacent and next-nearest pairs in a uniform waveguide bank."""
    _require_index(waveguides, "waveguides")
    if waveguides < 1:
        raise ValueError("waveguides must be >= 1")
    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,
    )

    if _rust_backend_enabled():
        analyzer = _backend_analyzer(
            "py_ph_analyze_crosstalk_bank", py_ph_analyze_crosstalk_bank
        )
        return dict(
            analyzer(
                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,
            )
        )

    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": total,
        "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)

Analyze per-pair crosstalk for arbitrary waveguide geometry.

Source code in src/sc_neurocore/optics/_photonic_crosstalk.py

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]:
    """Analyze per-pair crosstalk for arbitrary waveguide geometry."""
    pair_count = len(pair_indices)
    if len(gaps_nm) != pair_count or len(coupling_lengths_um) != pair_count:
        raise ValueError(
            f"pair_indices ({pair_count}), gaps_nm ({len(gaps_nm)}) and "
            f"coupling_lengths_um ({len(coupling_lengths_um)}) must be equal length"
        )
    for index, (pair_a, pair_b) in enumerate(pair_indices):
        _require_index(pair_a, f"pair_indices[{index}][0]")
        _require_index(pair_b, f"pair_indices[{index}][1]")
        if pair_a == pair_b:
            raise ValueError(f"pair_indices[{index}] must name two distinct waveguides")

    pairs = [
        WaveguidePair(
            gap_nm=gap,
            coupling_length_um=length,
            wavelength_nm=wavelength_nm,
            core_index=core_index,
            cladding_index=cladding_index,
        )
        for gap, length in zip(gaps_nm, coupling_lengths_um)
    ]
    if _rust_backend_enabled() and pair_count > 0:
        analyzer = _backend_analyzer(
            "py_ph_analyze_crosstalk_pairs", py_ph_analyze_crosstalk_pairs
        )
        return dict(
            analyzer(
                pairs_a=[pair_a for pair_a, _ in pair_indices],
                pairs_b=[pair_b for _, pair_b in pair_indices],
                gaps_nm=list(gaps_nm),
                lengths_um=list(coupling_lengths_um),
                wavelength_nm=wavelength_nm,
                core_index=core_index,
                cladding_index=cladding_index,
            )
        )

    return {
        "pair_a": [pair_a for pair_a, _ in pair_indices],
        "pair_b": [pair_b for _, pair_b in pair_indices],
        "gap_nm": list(gaps_nm),
        "coupling_length_um": list(coupling_lengths_um),
        "coupling_coefficient_per_um": [pair.coupling_coefficient for pair in pairs],
        "coupling_ratio": [pair.coupling_ratio for pair in pairs],
        "isolation_db": [pair.isolation_db for pair in pairs],
        "num_pairs": pair_count,
        "backend": "python",
    }

WaveguidePair dataclass

Physical contract for one pair of adjacent optical waveguides.

Source code in src/sc_neurocore/optics/_photonic_crosstalk.py

@dataclass
class WaveguidePair:
    """Physical contract for one pair of adjacent optical waveguides."""

    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

    def __post_init__(self) -> None:
        """Validate the coupled-mode domain before numerical evaluation."""
        _require_positive(self.waveguide_width_nm, "waveguide_width_nm")
        _require_non_negative(self.gap_nm, "gap_nm")
        _require_non_negative(self.coupling_length_um, "coupling_length_um")
        _require_positive(self.core_index, "core_index")
        _require_positive(self.cladding_index, "cladding_index")
        if self.core_index <= self.cladding_index:
            raise ValueError("core_index must be greater than cladding_index")
        _require_positive(self.wavelength_nm, "wavelength_nm")

    @property
    def effective_index_diff(self) -> float:
        """Return the Marcatili-form even/odd effective-index difference."""
        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:
        """Return coupling coefficient κ per micrometre."""
        return math.pi * self.effective_index_diff / (self.wavelength_nm * 1e-3)

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

    @property
    def isolation_db(self) -> float:
        """Return pair isolation in decibels with a 300 dB numeric ceiling."""
        ratio = self.coupling_ratio
        if ratio < 1e-15:
            return 300.0
        return -10.0 * math.log10(max(ratio, 1e-30))

effective_index_diff property

Return the Marcatili-form even/odd effective-index difference.

coupling_coefficient property

Return coupling coefficient κ per micrometre.

coupling_ratio property

Return power coupling ratio at the end of the parallel run.

isolation_db property

Return pair isolation in decibels with a 300 dB numeric ceiling.

__post_init__()

Validate the coupled-mode domain before numerical evaluation.

Source code in src/sc_neurocore/optics/_photonic_crosstalk.py

def __post_init__(self) -> None:
    """Validate the coupled-mode domain before numerical evaluation."""
    _require_positive(self.waveguide_width_nm, "waveguide_width_nm")
    _require_non_negative(self.gap_nm, "gap_nm")
    _require_non_negative(self.coupling_length_um, "coupling_length_um")
    _require_positive(self.core_index, "core_index")
    _require_positive(self.cladding_index, "cladding_index")
    if self.core_index <= self.cladding_index:
        raise ValueError("core_index must be greater than cladding_index")
    _require_positive(self.wavelength_nm, "wavelength_nm")

PhotonicEmitter

Emit photonic netlists in dependency order for a selected PDK.

Source code in src/sc_neurocore/optics/_photonic_emitter.py

class PhotonicEmitter:
    """Emit photonic netlists in dependency order for a selected PDK."""

    def __init__(self, target_pdk: str = "generic_si_photonics"):
        if not isinstance(target_pdk, str) or not target_pdk.strip():
            raise ValueError("target_pdk must be a non-empty string")
        self.target_pdk = target_pdk

    def _topological_sort(self, nodes: List[Any]) -> List[Any]:
        """Return nodes in stable dependency order and reject malformed graphs."""
        node_ids = [node.id for node in nodes]
        if len(set(node_ids)) != len(node_ids):
            raise ValueError("photonic IR node identifiers must be unique")
        outputs = [node.output for node in nodes]
        if len(set(outputs)) != len(outputs):
            raise ValueError("photonic IR node outputs must be unique")

        in_degree = {node.id: 0 for node in nodes}
        node_map = {node.id: node for node in nodes}
        adjacency: dict[str, list[str]] = {node.id: [] for node in nodes}
        output_to_id = {node.output: node.id for node in nodes}

        for node in nodes:
            for input_name in node.inputs:
                if input_name in output_to_id:
                    adjacency[output_to_id[input_name]].append(node.id)
                    in_degree[node.id] += 1

        queue = deque(node_id for node_id, degree in in_degree.items() if degree == 0)
        sorted_nodes: List[Any] = []
        while queue:
            current = queue.popleft()
            sorted_nodes.append(node_map[current])
            for neighbour in adjacency[current]:
                in_degree[neighbour] -= 1
                if in_degree[neighbour] == 0:
                    queue.append(neighbour)
        if len(sorted_nodes) != len(nodes):
            raise ValueError("photonic IR graph contains a dependency cycle")
        return sorted_nodes

    def emit_lumerical_netlist(self, ir_graph: Any) -> str:
        """Emit a Lumerical-compatible photonic netlist from an IR graph."""
        if not hasattr(ir_graph, "nodes"):
            raise TypeError("ir_graph must expose a nodes collection")
        sorted_nodes = self._topological_sort(ir_graph.nodes)
        netlist = ["# SC-NeuroCore Photonic Design", f"# PDK: {self.target_pdk}", ""]

        for node in sorted_nodes:
            if node.type == "SC_AND":
                if len(node.inputs) < 2:
                    raise ValueError(f"SC_AND node {node.id!r} requires two inputs")
                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":
                if not node.inputs:
                    raise ValueError(f"LIF_MEMBRANE node {node.id!r} requires an input")
                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")

        return "\n".join(netlist)

emit_lumerical_netlist(ir_graph)

Emit a Lumerical-compatible photonic netlist from an IR graph.

Source code in src/sc_neurocore/optics/_photonic_emitter.py

def emit_lumerical_netlist(self, ir_graph: Any) -> str:
    """Emit a Lumerical-compatible photonic netlist from an IR graph."""
    if not hasattr(ir_graph, "nodes"):
        raise TypeError("ir_graph must expose a nodes collection")
    sorted_nodes = self._topological_sort(ir_graph.nodes)
    netlist = ["# SC-NeuroCore Photonic Design", f"# PDK: {self.target_pdk}", ""]

    for node in sorted_nodes:
        if node.type == "SC_AND":
            if len(node.inputs) < 2:
                raise ValueError(f"SC_AND node {node.id!r} requires two inputs")
            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":
            if not node.inputs:
                raise ValueError(f"LIF_MEMBRANE node {node.id!r} requires an input")
            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")

    return "\n".join(netlist)

FDTD2DSolver

Two-dimensional TE Yee-grid solver with split-field Berenger PML.

Source code in src/sc_neurocore/optics/_photonic_fdtd.py

class FDTD2DSolver:
    """Two-dimensional TE Yee-grid solver with split-field Berenger PML."""

    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,
    ):
        _require_count(nx, "nx", minimum=3)
        _require_count(ny, "ny", minimum=3)
        _require_count(pml_layers, "pml_layers", minimum=1)
        if pml_layers >= min(nx, ny):
            raise ValueError("pml_layers must be smaller than both grid dimensions")
        _require_positive(dx_um, "dx_um")
        _require_positive(dy_um, "dy_um")
        _require_positive(dt_factor, "dt_factor")
        if dt_factor > 1.0:
            raise ValueError("dt_factor must not exceed the two-dimensional CFL scale factor of 1")

        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[Any, Any] = np.zeros((nx, ny), dtype=np.float64)
        self.ezy: np.ndarray[Any, Any] = np.zeros((nx, ny), dtype=np.float64)
        self.ez: np.ndarray[Any, Any] = np.zeros((nx, ny), dtype=np.float64)
        self.hx: np.ndarray[Any, Any] = np.zeros((nx, ny), dtype=np.float64)
        self.hy: np.ndarray[Any, Any] = np.zeros((nx, ny), dtype=np.float64)
        self.n_map: npt.NDArray[np.float64] = np.ones((nx, ny), dtype=np.float64)
        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 electric-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 on the material map."""
        _require_count(y_center, "y_center")
        _require_count(width_cells, "width_cells", minimum=1)
        _require_count(x_start, "x_start")
        if not 0 <= y_center < self.ny:
            raise ValueError(f"y_center {y_center} is outside a {self.ny}-cell grid")
        if x_end is not None:
            _require_count(x_end, "x_end")
        _require_positive(refractive_index, "refractive_index")
        if refractive_index < 1.0:
            raise ValueError(f"Invalid refractive index: {refractive_index}. Must be >= 1.0.")

        effective_end = self.nx if x_end is None else x_end
        effective_start = max(0, min(self.nx, x_start))
        effective_end = max(0, min(self.nx, effective_end))
        if effective_start >= effective_end:
            raise ValueError("x_start must be smaller than x_end after grid clipping")
        y_lo = max(0, min(self.ny, y_center - width_cells // 2))
        y_hi = max(0, min(self.ny, y_lo + width_cells))
        self.n_map[effective_start:effective_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 two-dimensional Gaussian electric-field source."""
        _require_count(x, "x")
        _require_count(y, "y")
        _require_count(sigma_cells, "sigma_cells", minimum=1)
        _require_positive(wavelength_nm, "wavelength_nm")
        _require_non_negative(amplitude, "amplitude")
        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 the TE simulation by ``n_steps`` timesteps."""
        _require_count(n_steps, "n_steps")
        if not np.all(np.isfinite(self.n_map)) or np.any(self.n_map <= 0):
            raise ValueError("Refractive index must be finite and > 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):
            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, :]
            )
            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:
        """Return total squared electromagnetic field energy."""
        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:
        """Return the electric field at one in-bounds grid point."""
        _require_count(x, "x")
        _require_count(y, "y")
        if x >= self.nx or y >= self.ny:
            raise ValueError(f"field point ({x}, {y}) out of bounds [{self.nx}, {self.ny}]")
        return float(self.ez[x, y])

    def cross_section(self, x: int) -> np.ndarray[Any, Any]:
        """Return an independent electric-field cross-section at ``x``."""
        _require_count(x, "x")
        if x >= self.nx:
            raise ValueError(f"cross-section index {x} out of bounds for nx={self.nx}")
        return self.ez[x, :].copy()

    def snapshot(self) -> Tuple[np.ndarray[Any, Any], np.ndarray[Any, Any], np.ndarray[Any, Any]]:
        """Return independent 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 on the material map.

Source code in src/sc_neurocore/optics/_photonic_fdtd.py

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 on the material map."""
    _require_count(y_center, "y_center")
    _require_count(width_cells, "width_cells", minimum=1)
    _require_count(x_start, "x_start")
    if not 0 <= y_center < self.ny:
        raise ValueError(f"y_center {y_center} is outside a {self.ny}-cell grid")
    if x_end is not None:
        _require_count(x_end, "x_end")
    _require_positive(refractive_index, "refractive_index")
    if refractive_index < 1.0:
        raise ValueError(f"Invalid refractive index: {refractive_index}. Must be >= 1.0.")

    effective_end = self.nx if x_end is None else x_end
    effective_start = max(0, min(self.nx, x_start))
    effective_end = max(0, min(self.nx, effective_end))
    if effective_start >= effective_end:
        raise ValueError("x_start must be smaller than x_end after grid clipping")
    y_lo = max(0, min(self.ny, y_center - width_cells // 2))
    y_hi = max(0, min(self.ny, y_lo + width_cells))
    self.n_map[effective_start:effective_end, y_lo:y_hi] = refractive_index

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

Inject a two-dimensional Gaussian electric-field source.

Source code in src/sc_neurocore/optics/_photonic_fdtd.py

def inject_source(
    self,
    x: int,
    y: int,
    wavelength_nm: float = 1550.0,
    amplitude: float = 1.0,
    sigma_cells: int = 10,
) -> None:
    """Inject a two-dimensional Gaussian electric-field source."""
    _require_count(x, "x")
    _require_count(y, "y")
    _require_count(sigma_cells, "sigma_cells", minimum=1)
    _require_positive(wavelength_nm, "wavelength_nm")
    _require_non_negative(amplitude, "amplitude")
    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 the TE simulation by n_steps timesteps.

Source code in src/sc_neurocore/optics/_photonic_fdtd.py

def step(self, n_steps: int = 1) -> None:
    """Advance the TE simulation by ``n_steps`` timesteps."""
    _require_count(n_steps, "n_steps")
    if not np.all(np.isfinite(self.n_map)) or np.any(self.n_map <= 0):
        raise ValueError("Refractive index must be finite and > 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):
        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, :]
        )
        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()

Return total squared electromagnetic field energy.

Source code in src/sc_neurocore/optics/_photonic_fdtd.py

def field_energy(self) -> float:
    """Return total squared electromagnetic field energy."""
    return float(np.sum(self.ez**2) + np.sum(self.hx**2) + np.sum(self.hy**2))

field_at_point(x, y)

Return the electric field at one in-bounds grid point.

Source code in src/sc_neurocore/optics/_photonic_fdtd.py

def field_at_point(self, x: int, y: int) -> float:
    """Return the electric field at one in-bounds grid point."""
    _require_count(x, "x")
    _require_count(y, "y")
    if x >= self.nx or y >= self.ny:
        raise ValueError(f"field point ({x}, {y}) out of bounds [{self.nx}, {self.ny}]")
    return float(self.ez[x, y])

cross_section(x)

Return an independent electric-field cross-section at x.

Source code in src/sc_neurocore/optics/_photonic_fdtd.py

def cross_section(self, x: int) -> np.ndarray[Any, Any]:
    """Return an independent electric-field cross-section at ``x``."""
    _require_count(x, "x")
    if x >= self.nx:
        raise ValueError(f"cross-section index {x} out of bounds for nx={self.nx}")
    return self.ez[x, :].copy()

snapshot()

Return independent copies of all field components.

Source code in src/sc_neurocore/optics/_photonic_fdtd.py

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

FDTDSolver

One-dimensional Yee-grid solver for waveguide co-simulation.

The solver applies a quadratic-ramp multiplicative absorbing boundary at each end. It is a bounded reference implementation for pulse propagation, dispersion, and loss checks; use :class:FDTD2DSolver when split-field Berenger PML is required.

Source code in src/sc_neurocore/optics/_photonic_fdtd.py

class FDTDSolver:
    """One-dimensional Yee-grid solver for waveguide co-simulation.

    The solver applies a quadratic-ramp multiplicative absorbing boundary at
    each end. It is a bounded reference implementation for pulse propagation,
    dispersion, and loss checks; use :class:`FDTD2DSolver` when split-field
    Berenger PML is required.
    """

    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,
    ):
        _require_count(grid_size, "grid_size", minimum=3)
        _require_count(boundary_cells, "boundary_cells", minimum=1)
        if boundary_cells > grid_size:
            raise ValueError("boundary_cells cannot exceed grid_size")
        _require_positive(dx_um, "dx_um")
        _require_positive(dt_factor, "dt_factor")
        if dt_factor > 1.0:
            raise ValueError("dt_factor must not exceed the one-dimensional CFL limit of 1")
        _require_positive(refractive_index, "refractive_index")

        self.grid_size = grid_size
        self.dx = dx_um * 1e-6
        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 non-negative propagation loss in decibels per centimetre."""
        _require_non_negative(loss_db_per_cm, "loss_db_per_cm")
        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."""
        _require_count(position, "position")
        if position >= self.grid_size:
            raise ValueError(f"position {position} is outside a {self.grid_size}-cell grid")
        _require_positive(wavelength_nm, "wavelength_nm")
        _require_non_negative(amplitude, "amplitude")
        _require_finite(phase, "phase")

        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."""
        _require_count(n_steps, "n_steps")
        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:
        """Return total squared electromagnetic field energy."""
        return float(np.sum(self.ez**2) + np.sum(self.hy**2))

    def snapshot(self) -> Tuple[np.ndarray[Any, Any], np.ndarray[Any, Any]]:
        """Return independent copies of the electric and magnetic fields."""
        return self.ez.copy(), self.hy.copy()

set_loss(loss_db_per_cm)

Set non-negative propagation loss in decibels per centimetre.

Source code in src/sc_neurocore/optics/_photonic_fdtd.py

def set_loss(self, loss_db_per_cm: float) -> None:
    """Set non-negative propagation loss in decibels per centimetre."""
    _require_non_negative(loss_db_per_cm, "loss_db_per_cm")
    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_fdtd.py

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."""
    _require_count(position, "position")
    if position >= self.grid_size:
        raise ValueError(f"position {position} is outside a {self.grid_size}-cell grid")
    _require_positive(wavelength_nm, "wavelength_nm")
    _require_non_negative(amplitude, "amplitude")
    _require_finite(phase, "phase")

    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_fdtd.py

def step(self, n_steps: int = 1) -> None:
    """Advance the simulation by ``n_steps`` timesteps."""
    _require_count(n_steps, "n_steps")
    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()

Return total squared electromagnetic field energy.

Source code in src/sc_neurocore/optics/_photonic_fdtd.py

def field_energy(self) -> float:
    """Return total squared electromagnetic field energy."""
    return float(np.sum(self.ez**2) + np.sum(self.hy**2))

snapshot()

Return independent copies of the electric and magnetic fields.

Source code in src/sc_neurocore/optics/_photonic_fdtd.py

def snapshot(self) -> Tuple[np.ndarray[Any, Any], np.ndarray[Any, Any]]:
    """Return independent copies of the electric and magnetic fields."""
    return self.ez.copy(), self.hy.copy()

MeepAdapter

Build and execute Meep waveguide simulations when Meep is installed.

Source code in src/sc_neurocore/optics/_photonic_meep.py

class MeepAdapter:
    """Build and execute Meep waveguide simulations when Meep is installed."""

    @staticmethod
    def is_available() -> bool:
        """Return whether the optional Meep dependency is importable."""
        try:
            import meep as meep_module

            return meep_module is not None
        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 serialisable Meep waveguide-geometry description."""
        if not isinstance(target, PhotonicTarget):
            raise TypeError("target must be a PhotonicTarget")
        _require_positive(waveguide_width_um, "waveguide_width_um")
        _require_positive(length_um, "length_um")
        _require_positive(substrate_index, "substrate_index")
        if substrate_index < 1.0:
            raise ValueError("substrate_index must be at least one")

        core_index = 3.48 if target.wavelength_nm > 1000 else 2.0
        wavelength_um = target.wavelength_nm / 1000.0
        frequency = 1.0 / wavelength_um
        return {
            "cell_size": [length_um, 3.0 * waveguide_width_um, 0],
            "resolution": 20,
            "sources": [
                {
                    "type": "ContinuousSource"
                    if target.modulation == OpticalModulation.PHASE
                    else "GaussianSource",
                    "frequency": frequency,
                    "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]:
        """Execute a real Meep simulation and return its transmission record.

        The method fails closed when Meep is unavailable. Earlier versions
        returned invented transmission values on that path, which could be
        mistaken for simulation evidence.
        """
        _require_positive(run_time, "run_time")
        if not isinstance(geometry, dict):
            raise TypeError("geometry must be a dictionary")
        _validate_geometry(geometry)
        if not MeepAdapter.is_available():
            raise ImportError(
                "Meep is not installed; install Meep before requesting a photonic simulation"
            )

        import meep as mp

        cell_size = geometry["cell_size"]
        resolution = geometry["resolution"]
        source_spec = geometry["sources"][0]
        geometry_spec = geometry["geometry"][0]
        source_factory = (
            mp.ContinuousSource if source_spec["type"] == "ContinuousSource" else mp.GaussianSource
        )

        frequency = source_spec["frequency"]
        sources = [
            mp.Source(
                source_factory(frequency=frequency),
                component=mp.Ez,
                center=mp.Vector3(*source_spec["center"]),
                size=mp.Vector3(*source_spec["size"]),
            )
        ]
        material = mp.Medium(index=geometry_spec["material_index"])
        geometry_objects = [
            mp.Block(
                size=mp.Vector3(geometry_spec["size"][0], geometry_spec["size"][1]),
                center=mp.Vector3(*geometry_spec["center"]),
                material=material,
            )
        ]
        simulation = mp.Simulation(
            cell_size=mp.Vector3(*cell_size),
            resolution=resolution,
            sources=sources,
            geometry=geometry_objects,
            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]),
        )
        transmission = simulation.add_flux(frequency, 0, 1, flux_region)
        simulation.run(until=run_time)
        flux_data = mp.get_fluxes(transmission)
        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

Return whether the optional Meep dependency is importable.

Source code in src/sc_neurocore/optics/_photonic_meep.py

@staticmethod
def is_available() -> bool:
    """Return whether the optional Meep dependency is importable."""
    try:
        import meep as meep_module

        return meep_module is not None
    except ImportError:
        return False

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

Build a serialisable Meep waveguide-geometry description.

Source code in src/sc_neurocore/optics/_photonic_meep.py

@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 serialisable Meep waveguide-geometry description."""
    if not isinstance(target, PhotonicTarget):
        raise TypeError("target must be a PhotonicTarget")
    _require_positive(waveguide_width_um, "waveguide_width_um")
    _require_positive(length_um, "length_um")
    _require_positive(substrate_index, "substrate_index")
    if substrate_index < 1.0:
        raise ValueError("substrate_index must be at least one")

    core_index = 3.48 if target.wavelength_nm > 1000 else 2.0
    wavelength_um = target.wavelength_nm / 1000.0
    frequency = 1.0 / wavelength_um
    return {
        "cell_size": [length_um, 3.0 * waveguide_width_um, 0],
        "resolution": 20,
        "sources": [
            {
                "type": "ContinuousSource"
                if target.modulation == OpticalModulation.PHASE
                else "GaussianSource",
                "frequency": frequency,
                "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

Execute a real Meep simulation and return its transmission record.

The method fails closed when Meep is unavailable. Earlier versions returned invented transmission values on that path, which could be mistaken for simulation evidence.

Source code in src/sc_neurocore/optics/_photonic_meep.py

@staticmethod
def run_simulation(geometry: Dict[str, Any], run_time: float = 50.0) -> Dict[str, Any]:
    """Execute a real Meep simulation and return its transmission record.

    The method fails closed when Meep is unavailable. Earlier versions
    returned invented transmission values on that path, which could be
    mistaken for simulation evidence.
    """
    _require_positive(run_time, "run_time")
    if not isinstance(geometry, dict):
        raise TypeError("geometry must be a dictionary")
    _validate_geometry(geometry)
    if not MeepAdapter.is_available():
        raise ImportError(
            "Meep is not installed; install Meep before requesting a photonic simulation"
        )

    import meep as mp

    cell_size = geometry["cell_size"]
    resolution = geometry["resolution"]
    source_spec = geometry["sources"][0]
    geometry_spec = geometry["geometry"][0]
    source_factory = (
        mp.ContinuousSource if source_spec["type"] == "ContinuousSource" else mp.GaussianSource
    )

    frequency = source_spec["frequency"]
    sources = [
        mp.Source(
            source_factory(frequency=frequency),
            component=mp.Ez,
            center=mp.Vector3(*source_spec["center"]),
            size=mp.Vector3(*source_spec["size"]),
        )
    ]
    material = mp.Medium(index=geometry_spec["material_index"])
    geometry_objects = [
        mp.Block(
            size=mp.Vector3(geometry_spec["size"][0], geometry_spec["size"][1]),
            center=mp.Vector3(*geometry_spec["center"]),
            material=material,
        )
    ]
    simulation = mp.Simulation(
        cell_size=mp.Vector3(*cell_size),
        resolution=resolution,
        sources=sources,
        geometry=geometry_objects,
        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]),
    )
    transmission = simulation.add_flux(frequency, 0, 1, flux_region)
    simulation.run(until=run_time)
    flux_data = mp.get_fluxes(transmission)
    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),
    }

OpticalModulation

Bases: Enum

Optical modulation scheme.

Source code in src/sc_neurocore/optics/_photonic_types.py

class OpticalModulation(Enum):
    """Optical modulation scheme."""

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

OpticalPulse dataclass

Single optical pulse with phase and amplitude.

Source code in src/sc_neurocore/optics/_photonic_types.py

@dataclass
class OpticalPulse:
    """Single optical pulse with phase and amplitude."""

    phase: float
    amplitude: float
    wavelength_nm: float
    duration_ps: float

    def __post_init__(self) -> None:
        """Validate the physical pulse boundary."""
        _require_finite(self.phase, "phase")
        _require_finite(self.amplitude, "amplitude")
        if not 0.0 <= self.amplitude <= 1.0:
            raise ValueError(f"amplitude must be in [0, 1], got {self.amplitude}")
        _require_positive(self.wavelength_nm, "wavelength_nm")
        _require_positive(self.duration_ps, "duration_ps")

__post_init__()

Validate the physical pulse boundary.

Source code in src/sc_neurocore/optics/_photonic_types.py

def __post_init__(self) -> None:
    """Validate the physical pulse boundary."""
    _require_finite(self.phase, "phase")
    _require_finite(self.amplitude, "amplitude")
    if not 0.0 <= self.amplitude <= 1.0:
        raise ValueError(f"amplitude must be in [0, 1], got {self.amplitude}")
    _require_positive(self.wavelength_nm, "wavelength_nm")
    _require_positive(self.duration_ps, "duration_ps")

PhotonicTarget dataclass

Hardware target specification for a photonic backend.

Source code in src/sc_neurocore/optics/_photonic_types.py

@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

    def __post_init__(self) -> None:
        """Validate target metadata before it reaches a compiler backend."""
        if not isinstance(self.name, str) or not self.name.strip():
            raise ValueError("name must be a non-empty string")
        if not isinstance(self.modulator_type, str) or not self.modulator_type.strip():
            raise ValueError("modulator_type must be a non-empty string")
        if not isinstance(self.modulation, OpticalModulation):
            raise TypeError("modulation must be an OpticalModulation")
        _require_positive(self.wavelength_nm, "wavelength_nm")
        _require_positive(self.q_factor, "q_factor")
        _require_non_negative(self.insertion_loss_db, "insertion_loss_db")
        _require_finite(self.thermo_optic_coeff, "thermo_optic_coeff")

    @classmethod
    def lightmatter(cls) -> PhotonicTarget:
        """Return a Lightmatter-style photonic target profile."""
        return cls("Lightmatter", 1550.0, OpticalModulation.PHASE, "MZI", 20000.0, 0.3)

    @classmethod
    def silicon_photonics(cls) -> PhotonicTarget:
        """Return a generic silicon-photonics target profile."""
        return cls("SiPh-Generic", 1310.0, OpticalModulation.AMPLITUDE, "Microring", 12000.0, 0.8)

    @classmethod
    def two_d_waveguide(cls) -> PhotonicTarget:
        """Return a two-dimensional-material waveguide target profile."""
        return cls("2D-Material", 850.0, OpticalModulation.HYBRID, "MZI", 5000.0, 1.2)

__post_init__()

Validate target metadata before it reaches a compiler backend.

Source code in src/sc_neurocore/optics/_photonic_types.py

def __post_init__(self) -> None:
    """Validate target metadata before it reaches a compiler backend."""
    if not isinstance(self.name, str) or not self.name.strip():
        raise ValueError("name must be a non-empty string")
    if not isinstance(self.modulator_type, str) or not self.modulator_type.strip():
        raise ValueError("modulator_type must be a non-empty string")
    if not isinstance(self.modulation, OpticalModulation):
        raise TypeError("modulation must be an OpticalModulation")
    _require_positive(self.wavelength_nm, "wavelength_nm")
    _require_positive(self.q_factor, "q_factor")
    _require_non_negative(self.insertion_loss_db, "insertion_loss_db")
    _require_finite(self.thermo_optic_coeff, "thermo_optic_coeff")

lightmatter() classmethod

Return a Lightmatter-style photonic target profile.

Source code in src/sc_neurocore/optics/_photonic_types.py

@classmethod
def lightmatter(cls) -> PhotonicTarget:
    """Return a Lightmatter-style photonic target profile."""
    return cls("Lightmatter", 1550.0, OpticalModulation.PHASE, "MZI", 20000.0, 0.3)

silicon_photonics() classmethod

Return a generic silicon-photonics target profile.

Source code in src/sc_neurocore/optics/_photonic_types.py

@classmethod
def silicon_photonics(cls) -> PhotonicTarget:
    """Return a generic silicon-photonics target profile."""
    return cls("SiPh-Generic", 1310.0, OpticalModulation.AMPLITUDE, "Microring", 12000.0, 0.8)

two_d_waveguide() classmethod

Return a two-dimensional-material waveguide target profile.

Source code in src/sc_neurocore/optics/_photonic_types.py

@classmethod
def two_d_waveguide(cls) -> PhotonicTarget:
    """Return a two-dimensional-material waveguide target profile."""
    return cls("2D-Material", 850.0, OpticalModulation.HYBRID, "MZI", 5000.0, 1.2)