Time-Delayed Coupling — Axonal Delay Kuramoto¶

The delay module implements Kuramoto dynamics with time-delayed coupling: oscillator \(j\)'s influence on oscillator \(i\) arrives with a discrete time delay \(\tau\). This models axonal propagation delays in neural circuits, signal latency in communication networks, and transport delays in distributed control systems.

Time-delayed coupling qualitatively changes the dynamics: it can stabilise otherwise unstable states, induce multistability, and generate travelling wave patterns not possible with instantaneous coupling.

1. Mathematical Formalism¶

1.1 The Delayed Kuramoto Equation¶

\[ \frac{d\theta_i}{dt} = \omega_i + \sum_{j=1}^{N} K_{ij} \sin(\theta_j(t - \tau) - \theta_i(t) - \alpha_{ij}) + \zeta \sin(\Psi - \theta_i(t)) \]

The key difference from standard Kuramoto: \(\theta_j\) in the coupling term is evaluated at time \(t - \tau\), not at the current time \(t\). This makes the system a delay differential equation (DDE), requiring the full history \(\theta(s)\) for \(s \in [t - \tau, t]\) as initial condition.

1.2 Discrete Delay Implementation¶

The module discretises the delay as an integer number of timesteps:

\[ \tau = \text{delay\_steps} \times \Delta t \]

Phase history is stored in a circular buffer (Python deque). At each step:

Current phases \(\theta(t)\) are appended to the buffer
Delayed phases \(\theta(t - \tau)\) are retrieved from \(\text{delay\_steps}\) positions back in the buffer
Coupling is computed using \(\theta_j^{\text{delayed}}\) and \(\theta_i^{\text{current}}\)
Euler integration advances \(\theta_i\)

\[ \theta_i^{n+1} = \left(\theta_i^n + \Delta t \left( \omega_i + \sum_j K_{ij} \sin(\theta_j^{n - d} - \theta_i^n - \alpha_{ij}) + \zeta\sin(\Psi - \theta_i^n) \right)\right) \bmod 2\pi \]

where \(d = \text{delay\_steps}\).

1.3 Buffer Initialisation¶

During the first \(d\) steps, insufficient history exists. The implementation falls back to current phases (\(\theta_j(t)\) instead of \(\theta_j(t-\tau)\)) until the buffer is filled. This is equivalent to assuming constant initial history \(\theta(s) = \theta(0)\) for \(s \in [-\tau, 0]\).

1.4 Delay Buffer (Standalone)¶

The DelayBuffer class provides a reusable circular buffer for phase history:

push(phases) — store a snapshot
get_delayed(delay_steps) — retrieve from \(d\) steps ago
clear() — reset buffer
length — current buffer occupancy

This can be used independently for custom delay schemes (variable delay, per-oscillator delay, etc.).

1.5 Effect of Delay on Synchronisation¶

For the standard Kuramoto model with uniform delay \(\tau\) and all-to-all coupling, the synchronised state satisfies (Yeung & Strogatz 1999):

\[ R = R_0 \cos(\omega_0 \tau) \]

where \(R_0\) is the undelayed order parameter and \(\omega_0\) is the mean frequency. Key consequences:

Constructive delay: when \(\omega_0 \tau = 2k\pi\) (integer multiples of the oscillation period), delay has no effect on \(R\)
Destructive delay: when \(\omega_0 \tau = (2k+1)\pi\), delay maximally reduces \(R\) and can destroy synchronisation
Critical delay: \(\tau_c = \pi / (2\omega_0)\) is the quarter-period threshold where \(R\) drops to zero

For non-uniform frequencies, the picture is more complex — delays can create chimera states (partial synchronisation with spatial structure).

1.6 Multistability from Delay¶

Time delay induces multistability: multiple stable phase-locked solutions coexist at the same coupling strength \(K\). The number of coexisting attractors grows with \(\tau\):

\[ N_{\text{attractors}} \sim \lfloor 2\omega_0 \tau / \pi \rfloor + 1 \]

This is a fundamental difference from instantaneous coupling, where (for identical oscillators) the synchronised state is unique.

1.7 Travelling Waves¶

For ring topologies with delay, travelling wave solutions exist:

\[ \theta_j(t) = \Omega t + j \cdot \phi_0 \]

where \(\Omega\) is the collective frequency (different from \(\omega_0\)) and \(\phi_0\) is the spatial wave number. The delay selects specific wave numbers — longer delays favour shorter wavelengths.

These travelling waves are observable in neural cortex (beta waves propagating across cortical sheets) and in SCPN Layer 4 (synchronisation dynamics).

1.8 Stability of Delayed Sync¶

The linearised stability of the synchronised state depends on the eigenvalues of the delayed coupling Jacobian:

\[ \lambda = -K R \cos(\alpha) e^{-\lambda \tau} + i\omega_0 \]

This transcendental equation (characteristic equation of the DDE) has infinitely many roots, unlike the ODE case. The dominant root determines stability. For small \(\tau\): \(\text{Re}(\lambda) \approx -K R \cos(\alpha)\) (same as undelayed). For large \(\tau\): \(\text{Re}(\lambda)\) oscillates, creating stability windows (Arnold tongues in the \((\tau, K)\) plane).

1.9 Delay Distribution¶

In biological networks, delays are not uniform — they follow a distribution \(p(\tau)\) depending on axon length, myelination, etc. The effective coupling with distributed delay is:

\[ \sum_j K_{ij} \int_0^\infty p(\tau) \sin(\theta_j(t-\tau) - \theta_i(t)) \, d\tau \]

SPO currently supports only fixed discrete delay. For distributed delays, sample multiple delay_steps values and weight the coupling accordingly.

2. Theoretical Context¶

2.1 Historical Background¶

Delayed coupling in oscillator networks was first studied systematically by Schuster & Wagner (1989) and Niebur, Schuster & Kammen (1991). Yeung & Strogatz (1999) provided the definitive analytical treatment of the delayed Kuramoto model, deriving the modified self-consistency equation and stability conditions.

The topic gained renewed interest with Dhamala, Jirsa & Ding (2004), who showed that time delays can enhance synchronisation in neural networks — contrary to naive expectation.

2.2 Applications¶

Domain	Source of delay	Typical \(\tau\)
Neural circuits	Axonal conduction velocity	1–50 ms
Brain regions	White matter tracts	5–200 ms
Power grids	Communication latency	10–500 ms
Distributed computing	Network RTT	1–100 ms
SCPN Layer 12	Gaian mesh UDP	10–100 ms

In SCPN, each layer has a characteristic timescale \(\tau_\ell\) (from SCPN_LAYER_TIMESCALES). Inter-layer coupling naturally has delays proportional to the timescale mismatch between layers.

2.3 Role in SPO Pipeline¶

The DelayedEngine extends the base UPDEEngine pattern for systems where instantaneous coupling is physically unrealistic. Common use:

Neural circuit simulations with axonal delays
Distributed SPO instances coupled via GaianMeshBridge (Layer 12)
Frequency-dependent delay modelling (higher frequencies → shorter effective delay)

2.4 Comparison with Other Delay Methods¶

Method	Delay type	Accuracy	Memory
`DelayedEngine` (Euler + buffer)	Fixed discrete	\(O(\Delta t)\)	\(O(d \cdot N)\)
DDE solver (e.g., `jitcdde`)	Continuous, interpolated	\(O(\Delta t^4)\)	\(O(d \cdot N)\)
State augmentation	Converts DDE to large ODE	\(O(\Delta t^4)\)	\(O(d \cdot N)\)

SPO uses fixed discrete delay for performance. For high-accuracy research requiring continuous delay, interface with jitcdde via the adapters.

3. Pipeline Position¶

┌──────────────┐     ┌──────────────────┐     ┌──────────────┐
│ coupling/    │────→│ DelayedEngine    │────→│ monitor/     │
│ knm.py       │     │                  │     │ order_params │
│ (K_ij)       │     │ Circular buffer  │     │ chimera      │
└──────────────┘     │ stores θ history │     └──────────────┘
                     │                  │
┌──────────────┐     │ θ_j(t-τ) used   │
│ oscillators/ │────→│ in coupling      │
│ base.py      │     └──────────────────┘
│ (ω_i)        │
└──────────────┘     ┌──────────────────┐
                     │ DelayBuffer      │
┌──────────────┐     │ (standalone)     │
│ adapters/    │────→│ push/get_delayed │
│ GaianMesh    │     └──────────────────┘
│ (network τ)  │
└──────────────┘

Inputs: - phases (N,) — current phases - omegas (N,) — natural frequencies - knm (N,N) — coupling matrix - alpha (N,N) — phase-lag (optional) - zeta, psi — external drive - delay_steps — discrete delay \(d\) (constructor parameter)

Outputs: - phases (N,) — updated phases after Euler step with delayed coupling

4. Features¶

4.1 Two Classes¶

Class	Purpose
`DelayBuffer`	Standalone circular buffer for phase history
`DelayedEngine`	Full integration engine with built-in delay buffer

4.2 Polyglot Acceleration¶

DelayedEngine.run() dispatches through a fastest-first backend chain resolved once at import (rust → mojo → julia → go → python). Each backend integrates the identical discrete delayed-Kuramoto recurrence over a ring buffer of delay_steps + 1 phase snapshots and returns the final phases wrapped into \([0, 2\pi)\):

Backend	Mechanism	Parity tolerance vs Python
`rust`	`spo-kernel` `delayed_kuramoto_run_rust` (PyO3)	\(10^{-9}\)
`mojo`	`mojo/delay_mojo` (stdin/stdout text protocol)	\(10^{-6}\)
`julia`	`julia/delay.jl` (`DelayJL.delayed_kuramoto_run`)	\(10^{-9}\)
`go`	`go/libdelay.so` `DelayedKuramotoRun` (cgo `c-shared`)	\(10^{-9}\)
`python`	`_python_run` NumPy ring-buffer reference	exact

ACTIVE_BACKEND names the resolved fastest backend; AVAILABLE_BACKENDS lists every backend that loaded successfully (Python is always present as the final floor). Each non-Python adapter validates its state arrays, scalars, and step counts via validate_delay_backend_inputs before any runtime loads, so malformed inputs fail closed without touching a foreign runtime. If a faster backend is missing the engine silently falls through to the next; results are bit-identical to the NumPy reference within the tolerances above.

Phase, frequency, coupling, and phase-lag arrays must contain numeric values, not stringified numbers. This check runs before NumPy float coercion on the public engine and direct Go/Julia/Mojo adapter boundaries. Shared backend outputs and Julia raw returns use the same rule before finite phase-vector, cardinality, and [0, 2*pi) checks. Mojo's stdout is intentionally textual; its protocol parser converts the declared output lines before applying the shared numeric output contract.

Note: DelayedEngine.step() is Python-only — the polyglot backends delegate at the run() level for batch efficiency, since the ring buffer must persist across the full step sequence.

4.3 Automatic Fallback¶

During initial \(d\) steps, the engine uses current phases instead of delayed phases. No special initialisation required — just call run().

4.4 Configurable Delay¶

delay_steps is set at construction time. The delay in physical units: \(\tau = d \times \Delta t\). For dt=0.01 and delay_steps=10: \(\tau = 0.1\) time units.

5. Usage Examples¶

5.1 Basic Delayed Coupling¶

import numpy as np
from scpn_phase_orchestrator.upde.delay import DelayedEngine
from scpn_phase_orchestrator.upde.order_params import compute_order_parameter

N = 32
rng = np.random.default_rng(42)
phases = rng.uniform(0, 2 * np.pi, N)
omegas = np.ones(N)  # identical frequencies
knm = np.full((N, N), 2.0 / N)
np.fill_diagonal(knm, 0.0)

# 10-step delay → τ = 10 × 0.01 = 0.1 time units
engine = DelayedEngine(N, dt=0.01, delay_steps=10)
phases = engine.run(phases, omegas, knm, n_steps=500)

R, _ = compute_order_parameter(phases)
print(f"Delayed R = {R:.4f}")

5.2 Delay Effect on Synchronisation¶

# Compare R for different delays
for d in [0, 5, 10, 20, 50]:
    eng = DelayedEngine(N, dt=0.01, delay_steps=max(d, 1))
    p = rng.uniform(0, 2 * np.pi, N).copy()
    p = eng.run(p, omegas, knm, n_steps=1000)
    R, _ = compute_order_parameter(p)
    tau = d * 0.01
    print(f"τ={tau:.2f}: R={R:.4f}")
# Expected: R decreases then oscillates with increasing τ

5.3 Constructive vs Destructive Delay¶

# For identical oscillators at ω₀=1: constructive at τ=2πk/ω₀, destructive at τ=π(2k+1)/ω₀
omega_0 = 1.0
omegas = np.ones(N) * omega_0

# Constructive: τ = 2π (one full period)
d_constructive = int(2 * np.pi / (omega_0 * 0.01))  # ~628 steps
eng_c = DelayedEngine(N, dt=0.01, delay_steps=d_constructive)
p_c = eng_c.run(rng.uniform(0, 2*np.pi, N), omegas, knm, n_steps=2000)
R_c, _ = compute_order_parameter(p_c)

# Destructive: τ = π (half period)
d_destructive = int(np.pi / (omega_0 * 0.01))  # ~314 steps
eng_d = DelayedEngine(N, dt=0.01, delay_steps=d_destructive)
p_d = eng_d.run(rng.uniform(0, 2*np.pi, N), omegas, knm, n_steps=2000)
R_d, _ = compute_order_parameter(p_d)

print(f"Constructive (τ=2π): R={R_c:.4f}")
print(f"Destructive  (τ=π):  R={R_d:.4f}")

5.4 Standalone DelayBuffer¶

from scpn_phase_orchestrator.upde.delay import DelayBuffer

buf = DelayBuffer(n_oscillators=8, max_delay_steps=20)

# Fill buffer with phase snapshots
for t in range(30):
    phases_t = np.sin(np.arange(8) * 0.1 + t * 0.05)
    buf.push(phases_t)

# Retrieve from 5 steps ago
delayed = buf.get_delayed(5)
print(f"Buffer length: {buf.length}")  # 20 (capped at max)

5.4 Neural Circuit with Axonal Delay¶

# Model cortical columns with 5ms axonal delay
dt = 0.001  # 1ms timestep
tau_axon = 0.005  # 5ms
delay_steps = int(tau_axon / dt)  # 5

N = 16
engine = DelayedEngine(N, dt=dt, delay_steps=delay_steps)
omegas = np.ones(N) * 40.0 * 2 * np.pi  # 40 Hz gamma band

phases = rng.uniform(0, 2 * np.pi, N)
phases = engine.run(phases, omegas, knm, n_steps=10000)
R, _ = compute_order_parameter(phases)
print(f"Gamma-band coherence with 5ms delay: R={R:.4f}")

5.5 Sakaguchi + Delay¶

alpha = np.full((N, N), np.pi / 8)
np.fill_diagonal(alpha, 0.0)

engine = DelayedEngine(N, dt=0.01, delay_steps=15)
phases = engine.run(phases, omegas, knm, alpha=alpha, n_steps=1000)
R, _ = compute_order_parameter(phases)
print(f"Frustrated + delayed: R={R:.4f}")

5.8 Delay Sweep (Arnold Tongue)¶

# Map R across (K, τ) parameter space
import numpy as np

N = 32
omegas = np.ones(N)
K_vals = np.linspace(0.5, 5.0, 20)
d_vals = range(1, 100, 5)  # delay in steps

R_map = np.zeros((len(K_vals), len(d_vals)))
for i, K in enumerate(K_vals):
    knm_k = np.full((N, N), K / N)
    np.fill_diagonal(knm_k, 0.0)
    for j, d in enumerate(d_vals):
        eng = DelayedEngine(N, dt=0.01, delay_steps=d)
        p = rng.uniform(0, 2 * np.pi, N)
        p = eng.run(p, omegas, knm_k, n_steps=1000)
        R_map[i, j], _ = compute_order_parameter(p)

# R_map shows Arnold tongue structure in (K, τ) plane

5.9 External Drive with Delay¶

# Entrainment with delayed feedback
engine = DelayedEngine(N, dt=0.01, delay_steps=20)
phases = engine.run(
    phases, omegas, knm,
    zeta=0.3, psi=np.pi/2,  # drive towards π/2
    n_steps=2000,
)
R, psi_mean = compute_order_parameter(phases)
print(f"Delayed entrainment: R={R:.4f}, ψ={psi_mean:.4f}")

6. Technical Reference¶

6.1 Module API¶

delay ¶

Time-delayed Kuramoto buffer and engine with validated phase history.

DelayBuffer stores copied finite phase snapshots in a bounded deque, and DelayedEngine advances phases with delayed coupling, optional external forcing, and Rust acceleration when available. Constructors and step inputs reject non-positive dimensions, non-finite scalars, shape-mismatched arrays, and boolean or numeric-string aliases before integration so delayed history never aliases invalid caller state.

Classes¶

DelayBuffer ¶

DelayBuffer(n_oscillators: int, max_delay_steps: int)

Circular buffer storing phase history for delayed coupling.

Stores last max_delay_steps snapshots. Retrieves phases from delay_steps steps ago.

Source code in src/scpn_phase_orchestrator/upde/delay.py

def __init__(self, n_oscillators: int, max_delay_steps: int):
    self._n = _validate_positive_int(n_oscillators, name="n_oscillators")
    self._max = _validate_positive_int(max_delay_steps, name="max_delay_steps")
    self._buffer: deque[FloatArray] = deque(maxlen=self._max)

Attributes¶

length `property` ¶

length: int

Number of snapshots currently stored.

Returns¶

int Number of snapshots currently stored.

Methods:¶

push ¶

push(phases: FloatArray) -> None

Append a phase snapshot to the buffer.

Parameters¶

phases : FloatArray Oscillator phases in radians, shape (N,).

Source code in src/scpn_phase_orchestrator/upde/delay.py

def push(self, phases: FloatArray) -> None:
    """Append a phase snapshot to the buffer.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    """
    phases64 = _validate_state_array(phases, name="phases", shape=(self._n,))
    self._buffer.append(phases64.copy())

get_delayed ¶

get_delayed(delay_steps: int) -> FloatArray | None

Return phases from delay_steps ago, or None if not enough history.

Parameters¶

delay_steps : int Number of steps in the past to retrieve from the delay buffer.

Returns¶

FloatArray | None The phase snapshot from delay_steps ago, or None if history is short.

Source code in src/scpn_phase_orchestrator/upde/delay.py

def get_delayed(self, delay_steps: int) -> FloatArray | None:
    """Return phases from `delay_steps` ago, or None if not enough history.

    Parameters
    ----------
    delay_steps : int
        Number of steps in the past to retrieve from the delay buffer.

    Returns
    -------
    FloatArray | None
        The phase snapshot from ``delay_steps`` ago, or ``None`` if history is
        short.
    """
    delay = _validate_positive_int(delay_steps, name="delay_steps")
    if delay > len(self._buffer):
        return None
    return self._buffer[-delay]

clear ¶

clear() -> None

Discard all stored phase snapshots.

Source code in src/scpn_phase_orchestrator/upde/delay.py

def clear(self) -> None:
    """Discard all stored phase snapshots."""
    self._buffer.clear()

DelayedEngine ¶

DelayedEngine(
    n_oscillators: int, dt: float, delay_steps: int = 1
)

Kuramoto with time-delayed coupling.

dθ_i/dt = ω_i + Σ_j K_ij sin(θ_j(t-τ) - θ_i(t) - α_ij)

Source code in src/scpn_phase_orchestrator/upde/delay.py

def __init__(self, n_oscillators: int, dt: float, delay_steps: int = 1):
    self._n = _validate_positive_int(n_oscillators, name="n_oscillators")
    self._dt = _validate_positive_float(dt, name="dt")
    self._delay_steps = _validate_positive_int(delay_steps, name="delay_steps")
    self._buffer: deque[FloatArray] = deque(maxlen=self._delay_steps + 1)

Attributes¶

delay_steps `property` ¶

delay_steps: int

Return the configured discrete coupling delay.

Returns¶

int Return the configured discrete coupling delay.

Methods:¶

step ¶

step(
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float = 0.0,
    psi: float = 0.0,
    alpha: FloatArray | None = None,
    step_idx: int = 0,
) -> FloatArray

Advance one delayed Kuramoto timestep from validated state arrays.

Parameters¶

phases : FloatArray Oscillator phases in radians, shape (N,). omegas : FloatArray Natural frequencies in rad/s, shape (N,). knm : FloatArray Coupling matrix K_nm, shape (N, N). zeta : float External drive strength ζ. psi : float External drive reference phase Ψ in radians. alpha : FloatArray | None Phase-lag matrix in radians, shape (N, N), or None for no lag. step_idx : int Zero-based index of the current step, used to address delayed coupling history.

Returns¶

FloatArray The phases after one delayed Kuramoto step, in [0, 2π).

Source code in src/scpn_phase_orchestrator/upde/delay.py

def step(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float = 0.0,
    psi: float = 0.0,
    alpha: FloatArray | None = None,
    step_idx: int = 0,
) -> FloatArray:
    """Advance one delayed Kuramoto timestep from validated state arrays.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    zeta : float
        External drive strength ``ζ``.
    psi : float
        External drive reference phase ``Ψ`` in radians.
    alpha : FloatArray | None
        Phase-lag matrix in radians, shape ``(N, N)``, or ``None`` for no lag.
    step_idx : int
        Zero-based index of the current step, used to address delayed coupling
        history.

    Returns
    -------
    FloatArray
        The phases after one delayed Kuramoto step, in ``[0, 2π)``.
    """
    del step_idx
    phases64 = _validate_state_array(phases, name="phases", shape=(self._n,))
    omegas64 = _validate_state_array(omegas, name="omegas", shape=(self._n,))
    knm64 = _validate_state_array(knm, name="knm", shape=(self._n, self._n))
    alpha64: FloatArray
    if alpha is None:
        alpha64 = np.zeros((self._n, self._n), dtype=np.float64)
    else:
        alpha64 = _validate_state_array(
            alpha,
            name="alpha",
            shape=(self._n, self._n),
        )
    zeta = _validate_finite_float(zeta, name="zeta")
    psi = _validate_finite_float(psi, name="psi")
    self._buffer.append(phases64.copy())
    delayed = self._buffer[0] if len(self._buffer) > self._delay_steps else phases64
    diff = delayed[np.newaxis, :] - phases64[:, np.newaxis] - alpha64
    coupling = np.sum(knm64 * np.sin(diff), axis=1)
    dtheta = omegas64 + coupling
    if zeta != 0.0:
        dtheta += zeta * np.sin(psi - phases64)
    step_out: FloatArray = (phases64 + self._dt * dtheta) % TWO_PI
    return step_out

run ¶

run(
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float = 0.0,
    psi: float = 0.0,
    alpha: FloatArray | None = None,
    n_steps: int = 100,
) -> FloatArray

Run delayed Kuramoto integration for n_steps validated steps.

Parameters¶

phases : FloatArray Oscillator phases in radians, shape (N,). omegas : FloatArray Natural frequencies in rad/s, shape (N,). knm : FloatArray Coupling matrix K_nm, shape (N, N). zeta : float External drive strength ζ. psi : float External drive reference phase Ψ in radians. alpha : FloatArray | None Phase-lag matrix in radians, shape (N, N), or None for no lag. n_steps : int Number of integration steps to run.

Returns¶

FloatArray The final phases after n_steps delayed Kuramoto steps.

Source code in src/scpn_phase_orchestrator/upde/delay.py

def run(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    knm: FloatArray,
    zeta: float = 0.0,
    psi: float = 0.0,
    alpha: FloatArray | None = None,
    n_steps: int = 100,
) -> FloatArray:
    """Run delayed Kuramoto integration for ``n_steps`` validated steps.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    knm : FloatArray
        Coupling matrix ``K_nm``, shape ``(N, N)``.
    zeta : float
        External drive strength ``ζ``.
    psi : float
        External drive reference phase ``Ψ`` in radians.
    alpha : FloatArray | None
        Phase-lag matrix in radians, shape ``(N, N)``, or ``None`` for no lag.
    n_steps : int
        Number of integration steps to run.

    Returns
    -------
    FloatArray
        The final phases after ``n_steps`` delayed Kuramoto steps.
    """
    n_steps = _validate_positive_int(n_steps, name="n_steps")
    phases64 = _validate_state_array(phases, name="phases", shape=(self._n,))
    omegas64 = _validate_state_array(omegas, name="omegas", shape=(self._n,))
    knm64 = _validate_state_array(knm, name="knm", shape=(self._n, self._n))
    alpha64: FloatArray
    if alpha is None:
        alpha64 = np.zeros((self._n, self._n), dtype=np.float64)
    else:
        alpha64 = _validate_state_array(
            alpha,
            name="alpha",
            shape=(self._n, self._n),
        )
    zeta = _validate_finite_float(zeta, name="zeta")
    psi = _validate_finite_float(psi, name="psi")
    knm_flat = np.ascontiguousarray(knm64.ravel(), dtype=np.float64)
    alpha_flat = np.ascontiguousarray(alpha64.ravel(), dtype=np.float64)
    backend_fn = _dispatch()
    if backend_fn is not None:
        try:
            return _validate_phase_output(
                backend_fn(
                    phases64,
                    omegas64,
                    knm_flat,
                    alpha_flat,
                    self._n,
                    zeta,
                    psi,
                    self._dt,
                    self._delay_steps,
                    n_steps,
                ),
                n_oscillators=self._n,
            )
        except (ImportError, RuntimeError, OSError, KeyError, ValueError):
            pass
    return _python_run(
        phases64,
        omegas64,
        knm_flat,
        alpha_flat,
        self._n,
        zeta,
        psi,
        self._dt,
        self._delay_steps,
        n_steps,
    )

Functions:¶

6.2 DelayBuffer API¶

Method	Description
`__init__(n_oscillators, max_delay_steps)`	Create buffer
`push(phases)`	Append snapshot (copies input)
`get_delayed(delay_steps) → NDArray or None`	Retrieve from \(d\) steps back
`length → int`	Current occupancy
`clear()`	Reset buffer

6.3 DelayedEngine API¶

Method/Property	Description
`__init__(n_oscillators, dt, delay_steps=1)`	Create engine
`delay_steps → int`	Configured delay
`step(phases, omegas, knm, zeta, psi, alpha, step_idx)`	Single step (Python)
`run(phases, omegas, knm, zeta, psi, alpha, n_steps)`	Batch (Rust or Python)

6.4 Parameters¶

Parameter	Type	Default	Description
`n_oscillators`	`int`	—	Number of oscillators
`dt`	`float`	—	Timestep
`delay_steps`	`int`	`1`	Discrete delay \(d\) (\(\tau = d \cdot \Delta t\))

7. Performance Benchmarks¶

7.1 Delay Overhead¶

The delay buffer adds minimal overhead: one deque.append() and one deque[index] per step — both \(O(1)\) regardless of buffer size.

N	delay_steps	Python (µs/step)	Rust run() (µs/step)	Speedup
16	10	42	7.1	5.9x
64	10	340	52	6.5x
256	10	8200	1100	7.5x
64	50	340	52	6.5x

Delay steps do NOT affect per-step cost significantly — the buffer is \(O(1)\) access. The dominant cost remains the \(O(N^2)\) coupling sum.

7.2 Memory¶

Buffer memory: \(d \times N \times 8\) bytes (float64).

N	delay_steps	Buffer size
16	10	1.3 KB
64	50	25 KB
256	100	200 KB
1024	200	1.6 MB

7.3 Complexity¶

Operation	Time	Space
`step()`	\(O(N^2)\)	\(O(N^2)\) scratch + \(O(d \cdot N)\) buffer
`run(n_steps)`	\(O(n \cdot N^2)\)	\(O(d \cdot N)\) buffer

7.4 Delay vs Undelayed Comparison¶

For the same system (N=64, K=3.0/N, 1000 steps):

Engine	delay_steps	µs/step (Rust)	R_final
`UPDEEngine` (Euler)	0	57	0.908
`DelayedEngine`	1	58	0.905
`DelayedEngine`	10	59	0.891
`DelayedEngine`	50	60	0.742
`DelayedEngine`	100	61	0.651

Overhead of delay buffer: < 5% at N=64. Physics effect is significant: \(R\) drops from 0.91 to 0.65 with \(\tau = 1.0\) time unit.

7.5 Recommended Settings¶

Application	dt	delay_steps	τ (time units)
Neural gamma (40 Hz)	0.001	5	5 ms
Neural beta (20 Hz)	0.001	25	25 ms
SCPN inter-layer	0.01	10–50	0.1–0.5
Gaian mesh (UDP)	0.01	5–20	0.05–0.2
Power grid comms	0.01	50–500	0.5–5.0

7.6 Polyglot Parity Gate¶

benchmarks/delay_benchmark.py (benchmark_delay_polyglot_parity_gate) records every declared backend's final phases against the NumPy reference and asserts two contracts per backend: with zero coupling and zero forcing the integrator must reduce to pure rotation \((\theta + \omega\,\Delta t\,n_{\text{steps}}) \bmod 2\pi\), and all output phases must lie in \([0, 2\pi)\). The gate is registered in the aggregate reference suite under the key delay_polyglot.

Representative run (n=16, delay_steps=3, n_steps=80, calls=1, seed=2026; evidence kind local_regression_non_isolated, no production timing claim — ML350 Gen8, single warmed call):

Backend	max abs error vs Python	ms/call	parity
`rust`	4.4e-16	0.93	✓
`mojo`	4.4e-16	22.64	✓
`julia`	4.4e-16	0.68	✓
`go`	8.9e-16	0.54	✓
`python`	0.0 (reference)	1.13	✓

Pure-rotation contract error: 8.9e-15 (≤ 1e-9 threshold). Single-call timings are noisy and not isolated; they exist to confirm every backend runs, not to rank production throughput.

7.7 When to Use DelayedEngine vs UPDEEngine¶

Use DelayedEngine when: - Physical delay is ≥ 1 timestep (\(\tau \geq \Delta t\)) - Delay is a significant fraction of the oscillation period (\(\tau \geq T/10\) where \(T = 2\pi/\omega_0\)) - You need to model travelling waves or delay-induced multistability

Use UPDEEngine when: - Coupling is effectively instantaneous (\(\tau \ll \Delta t\)) - The system is fast enough that delay is negligible

8. Citations¶

Yeung M.K.S., Strogatz S.H. (1999). Time delay in the Kuramoto model of coupled oscillators. Physical Review Letters 82(3):648–651. doi:10.1103/PhysRevLett.82.648
Schuster H.G., Wagner P. (1989). Mutual entrainment of two limit cycle oscillators with time delayed coupling. Progress of Theoretical Physics 81(5):939–945. doi:10.1143/PTP.81.939
Niebur E., Schuster H.G., Kammen D.M. (1991). Collective frequencies and metastability in networks of limit-cycle oscillators with time delay. Physical Review Letters 67(20):2753–2756. doi:10.1103/PhysRevLett.67.2753
Dhamala M., Jirsa V.K., Ding M. (2004). Enhancement of neural synchrony by time delay. Physical Review Letters 92(7):074104. doi:10.1103/PhysRevLett.92.074104
Kuramoto Y., Battogtokh D. (2002). Coexistence of coherence and incoherence in nonlocally coupled phase oscillators. Nonlinear Phenomena in Complex Systems 5(4):380–385.

Test Coverage¶

tests/test_delay.py — 20 tests: DelayBuffer push/get/clear/length, DelayedEngine step() shape, run() convergence, delay effect on R, zero delay equivalence, Sakaguchi + delay, external drive, dispatch fallback parity
tests/test_delay_backends.py — 46 tests: per-backend final-phase parity (Rust/Julia/Go at 1e-9, Mojo at 1e-6), cross-backend consistency, direct adapter boundary contracts (dtype, finiteness, real-valued, length, range validation), and loader/dispatch resolution
tests/test_delay_benchmark.py — 3 tests: polyglot parity gate slot reporting, control validation, and per-backend wall-clock timings
tests/test_reference_benchmark_suite.py — delay_polyglot gate slot and aggregate-suite key membership

Source¶

Python: src/scpn_phase_orchestrator/upde/delay.py (dispatcher + _python_run NumPy reference)
Validation: src/scpn_phase_orchestrator/experimental/accelerators/upde/_delay_validation.py
Rust: spo-kernel/crates/spo-engine/src/delay.rs
FFI: spo-kernel/crates/spo-ffi/src/lib.rs (delayed_kuramoto_run_rust)
Julia: julia/delay.jl (DelayJL.delayed_kuramoto_run)
Go: go/delay.go (DelayedKuramotoRun, c-shared → libdelay.so)
Mojo: mojo/delay.mojo (stdin/stdout text protocol → delay_mojo)
Backend adapters: src/scpn_phase_orchestrator/experimental/accelerators/upde/_delay_{julia,go,mojo}.py
Benchmark: benchmarks/delay_benchmark.py

Time-Delayed Coupling — Axonal Delay Kuramoto¶

1. Mathematical Formalism¶

1.1 The Delayed Kuramoto Equation¶

1.2 Discrete Delay Implementation¶

1.3 Buffer Initialisation¶

1.4 Delay Buffer (Standalone)¶

1.5 Effect of Delay on Synchronisation¶

1.6 Multistability from Delay¶

1.7 Travelling Waves¶

1.8 Stability of Delayed Sync¶

1.9 Delay Distribution¶

2. Theoretical Context¶

2.1 Historical Background¶

2.2 Applications¶

2.3 Role in SPO Pipeline¶

2.4 Comparison with Other Delay Methods¶

3. Pipeline Position¶

4. Features¶

4.1 Two Classes¶

4.2 Polyglot Acceleration¶

4.3 Automatic Fallback¶

4.4 Configurable Delay¶

5. Usage Examples¶

5.1 Basic Delayed Coupling¶

5.2 Delay Effect on Synchronisation¶

5.3 Constructive vs Destructive Delay¶

5.4 Standalone DelayBuffer¶

5.4 Neural Circuit with Axonal Delay¶

5.5 Sakaguchi + Delay¶

5.8 Delay Sweep (Arnold Tongue)¶

5.9 External Drive with Delay¶

6. Technical Reference¶

6.1 Module API¶

delay ¶

Classes¶

DelayBuffer ¶

Attributes¶

length property ¶

Returns¶

Methods:¶

push ¶

Parameters¶

get_delayed ¶

Parameters¶

Returns¶

clear ¶

DelayedEngine ¶

Attributes¶

delay_steps property ¶

Returns¶

Methods:¶

step ¶

Parameters¶

Returns¶

run ¶

Parameters¶

Returns¶

Functions:¶

6.2 DelayBuffer API¶

6.3 DelayedEngine API¶

6.4 Parameters¶

7. Performance Benchmarks¶

7.1 Delay Overhead¶

7.2 Memory¶

7.3 Complexity¶

7.4 Delay vs Undelayed Comparison¶

7.5 Recommended Settings¶

7.6 Polyglot Parity Gate¶

7.7 When to Use DelayedEngine vs UPDEEngine¶

8. Citations¶

Test Coverage¶

Source¶

length `property` ¶

delay_steps `property` ¶