Ott-Antonsen Mean-Field Reduction¶

1. Mathematical Formalism¶

The Ott-Antonsen Ansatz¶

For an infinite population of globally coupled Kuramoto oscillators with Lorentzian (Cauchy) natural frequency distribution:

\[g(\omega) = \frac{\Delta / \pi}{(\omega - \omega_0)^2 + \Delta^2}\]

where $\omega_0$ is the centre frequency and $\Delta$ is the half-width at half-maximum, the exact low-dimensional dynamics of the Kuramoto order parameter $z = R e^{i\psi}$ are given by:

\[\frac{dz}{dt} = -(\Delta + i\omega_0)\,z + \frac{K}{2}\left(z - |z|^2 z\right)\]

This is the Ott-Antonsen reduction — a single complex ODE that captures the exact macroscopic dynamics of the infinite-$N$ Kuramoto model.

Decomposition of the ODE¶

Separating real and imaginary parts $z = x + iy$:

\[\dot{x} = -\Delta x + \omega_0 y + \frac{K}{2}(1 - x^2 - y^2)\,x\]

\[\dot{y} = -\Delta y - \omega_0 x + \frac{K}{2}(1 - x^2 - y^2)\,y\]

The three terms have clear physical meanings: 1. $-\Delta z$: Damping from frequency spread (wider distribution → faster decay) 2. $-i\omega_0 z$: Rotation at the mean frequency 3. $\frac{K}{2}(z - |z|^2 z)$: Cubic self-coupling from the Lorentzian closure (Ott & Antonsen 2008, Eq. 9)

Steady-State Solution¶

Setting $\dot{R} = 0$ in the polar representation:

\[R_{ss} = \begin{cases} 0 & \text{if } K \leq K_c = 2\Delta \\ \sqrt{1 - \frac{2\Delta}{K}} & \text{if } K > K_c \end{cases}\]

This is exact for the Lorentzian distribution. The transition at $K_c = 2\Delta$ is continuous (second-order), with the universal Kuramoto scaling $R \propto \sqrt{K - K_c}$ near the transition.

Critical Coupling¶

\[K_c = 2\Delta\]

This is the Dörfler-Bullo criterion for globally coupled oscillators with Lorentzian spread. Larger frequency spread (larger $\Delta$) requires stronger coupling to achieve synchronisation.

Integration Method¶

RK4 (Runge-Kutta 4th order) on the complex ODE:

\[z_{n+1} = z_n + \frac{\Delta t}{6}(k_1 + 2k_2 + 2k_3 + k_4)\]

where $k_i$ are the standard RK4 stages applied to $f(z) = -(\Delta + i\omega_0)z + \frac{K}{2}(z - |z|^2 z)$.

2. Theoretical Context¶

Why Mean-Field Reduction?¶

Simulating $N$ oscillators costs $O(N^2)$ per step (all-to-all coupling). For $N = 10{,}000$ with 1000 steps, that is $10^{11}$ operations. The Ott-Antonsen reduction replaces this with a single complex ODE — constant cost per step, independent of $N$.

This is not an approximation: for Lorentzian $g(\omega)$ and $N \to \infty$, the reduction is exact. For finite $N$, it provides a mean-field prediction that converges as $1/\sqrt{N}$.

Domain of Validity¶

The Ott-Antonsen ansatz applies when: 1. Coupling is global (all-to-all): $K_{ij} = K/N$ for all $i, j$ 2. Frequency distribution is Lorentzian (or can be well approximated) 3. No noise (deterministic dynamics) 4. No phase frustration ($\alpha_{ij} = 0$)

For non-Lorentzian distributions, the ansatz is no longer exact but remains a useful approximation. For Gaussian $g(\omega)$, the error is typically $< 5\%$ in $R_{ss}$ for moderate coupling.

Historical Context¶

Kuramoto (1975): Original model, heuristic self-consistency for $R$
Strogatz (2000): Mean-field theory for finite $N$
Ott & Antonsen (2008): Exact low-dimensional reduction for Lorentzian $g(\omega)$
Ott & Antonsen (2009): Extension to multivariate coupling
Marvel, Mirollo & Strogatz (2009): Identified the OA manifold as the Möbius group (Watanabe-Strogatz theory)
Pikovsky & Rosenblum (2008): Analogous reduction via circular cumulants

The Ott-Antonsen result is considered one of the most significant theoretical advances in coupled oscillator theory since Kuramoto's original work.

Applications in SPO¶

The OttAntonsenReduction serves as a fast diagnostic: - Before running expensive N-oscillator simulations, predict whether the given $(K, \Delta)$ will produce synchronisation - Estimate the steady-state $R$ for parameter sweeps in seconds rather than hours - Validate N-oscillator results against the exact mean-field solution - Provide initial conditions for the full simulation (seed phases near the predicted $R_{ss}$)

Lorentzian Fitting¶

The predict_from_oscillators() method automatically estimates $(\omega_0, \Delta)$ from a set of measured frequencies: - $\omega_0 = \text{median}(\omega)$ (robust central tendency) - $\Delta = \text{IQR}(\omega) / 2$ (the IQR of a Lorentzian equals $2\Delta$)

This enables one-call diagnostics: pass frequencies, get predicted $R$.

3. Pipeline Position¶

Oscillators.extract() ──→ ω₁, ω₂, ..., ωₙ
                                │
                                ↓
     ┌── OttAntonsenReduction(ω₀, Δ, K, dt) ──┐
     │                                          │
     │  Input:  z₀ (seed), n_steps             │
     │  Param:  ω₀, Δ (from g(ω)), K, dt      │
     │  Method: RK4 (Rust/Mojo/Julia/Go/Python)│
     │  Output: OAState(z, R, ψ, K_c)          │
     │                                          │
     │  OR: predict_from_oscillators(ω, K)      │
     │      → fit Lorentzian → run → OAState    │
     │                                          │
     └──────────────────────────────────────────┘
                                │
                                ↓
              R_predicted vs R_measured (from UPDEEngine)
              → validation / parameter guidance

Input Contracts¶

Parameter	Type	Range	Source
`omega_0`	`float`	any	Median of natural frequencies
`delta`	`float`	$\geq 0$	Half-width of Lorentzian fit
`K`	`float`	finite real	Coupling strength; negative values model repulsive coupling
`dt`	`float`	$> 0$	Integration time step (default 0.01)
`z0`	`complex`	$	z_0
`n_steps`	`int`	$\geq 1$	Number of RK4 steps
`omegas`	real `FloatArray`	finite, non-empty, one-dimensional	Empirical sample for `predict_from_oscillators()`; boolean, complex, and numeric-string aliases are rejected before fitting

Output Contract¶

Field	Type	Range	Meaning
`z`	`complex`	$	z
`R`	`float`	$[0, 1]$	$
`psi`	`float`	$(-\pi, \pi]$	$\arg(z)$, mean phase
`K_c`	`float`	$\geq 0$	Critical coupling $2\Delta$

Accelerator Output Contract¶

Direct Go, Julia, Mojo, and Rust-backed public dispatch outputs must preserve the same Ott-Antonsen state invariant before an OAState is published. The backend tuple must contain finite non-boolean real scalars, the complex state must remain inside the OA unit disk, R must equal abs(z), and psi must match atan2(Im(z), Re(z)) for non-zero radius. Loader or runtime unavailability still falls through to the Python floor, but physics-contract faults propagate as validation errors. The Rust kernel runner validates its scalar inputs before conversion, and the optional Rust steady-state helper must return a finite real scalar inside [0, 1] before the public result is emitted.

4. Features¶

Exact reduction for Lorentzian $g(\omega)$ — not an approximation
RK4 integration — 4th-order accurate on the complex ODE
Rust FFI acceleration — 38-96x speedup over Python loop
Mojo subprocess boundary — accepts exactly four finite scalar records from the compiled OA runner before constructing an OAState
Public backend output replay — optional backend outputs are checked for R == |z|, psi == atan2(Im(z), Re(z)), and non-boolean finite scalars before publication
Typed empirical fitting boundary — measured-frequency samples reject boolean, complex, and numeric-string aliases before Lorentzian fitting
Analytical steady-state — steady_state_R() returns exact $R_{ss}$
Automatic Lorentzian fitting — predict_from_oscillators() fits $(\omega_0, \Delta)$ from measured frequencies
Critical coupling — K_c property returns $2\Delta$
Single-step and batch — step() for interactive, run() for batch
Immutable state — returns OAState dataclass, engine is reusable

5. Usage Examples¶

Basic: Predict R for Given Parameters¶

from scpn_phase_orchestrator.upde.reduction import OttAntonsenReduction

oa = OttAntonsenReduction(omega_0=0.0, delta=0.5, K=3.0, dt=0.01)

print(f"K_c = {oa.K_c:.2f}")           # 1.0
print(f"R_ss = {oa.steady_state_R():.4f}")  # sqrt(1 - 2*0.5/3.0) = 0.8165

state = oa.run(z0=complex(0.01, 0.0), n_steps=2000)
print(f"R (numerical) = {state.R:.4f}")  # Should match R_ss

Fast Diagnostic from Oscillator Frequencies¶

import numpy as np
from scpn_phase_orchestrator.upde.reduction import OttAntonsenReduction

# Given N oscillator frequencies (e.g. from EEG channels)
omegas = np.random.default_rng(42).standard_cauchy(100) * 0.3 + 1.0

oa = OttAntonsenReduction(omega_0=0, delta=0.1, K=1.0, dt=0.01)
prediction = oa.predict_from_oscillators(omegas, K=2.0)

print(f"Predicted R = {prediction.R:.4f}")
print(f"K_c = {prediction.K_c:.4f}")
print(f"Supercritical: {2.0 > prediction.K_c}")

Parameter Sweep¶

import numpy as np
from scpn_phase_orchestrator.upde.reduction import OttAntonsenReduction

delta = 0.5
K_values = np.linspace(0, 5, 100)

for K in K_values:
    oa = OttAntonsenReduction(omega_0=0, delta=delta, K=K, dt=0.01)
    R_analytical = oa.steady_state_R()
    state = oa.run(complex(0.01, 0.0), 3000)
    print(f"K={K:.2f}: R_analytical={R_analytical:.4f}, R_numerical={state.R:.4f}")

Validation Against Full Simulation¶

import numpy as np
from scpn_phase_orchestrator.upde.reduction import OttAntonsenReduction
from scpn_phase_orchestrator.upde.engine import UPDEEngine
from scpn_phase_orchestrator.upde.order_params import compute_order_parameter

N = 1000
K = 3.0
delta = 0.5

# OA prediction (milliseconds)
oa = OttAntonsenReduction(omega_0=0, delta=delta, K=K, dt=0.01)
R_oa = oa.steady_state_R()

# Full simulation (seconds)
omegas = np.random.default_rng(42).standard_cauchy(N) * delta
knm = np.full((N, N), K / N); np.fill_diagonal(knm, 0.0)
phases = np.random.default_rng(42).uniform(0, 2 * np.pi, N)

engine = UPDEEngine(N, dt=0.01)
for _ in range(5000):
    phases = engine.step(phases, omegas, knm, 0.0, 0.0, np.zeros((N, N)))
R_sim, _ = compute_order_parameter(phases)

print(f"OA prediction: R = {R_oa:.4f}")
print(f"Full sim (N={N}): R = {R_sim:.4f}")
# Error: typically < 0.05 for N >= 100

6. Technical Reference¶

Class: OttAntonsenReduction¶

reduction ¶

Exact mean-field (Ott-Antonsen) reduction for globally-coupled Kuramoto.

Uses a Lorentzian g(ω) and a 5-backend fallback chain.

Dynamics¶

On the Ott-Antonsen manifold the full N-oscillator Kuramoto system reduces to a single complex-scalar ODE:

dz/dt = −(Δ + iω₀)·z + (K/2)·(z − |z|²·z)

with z = R·e^{iψ} the mean-field order parameter, Δ the half-width of the Lorentzian g(ω), ω₀ its centre, and K the coupling strength.

Steady-state R_ss = √(1 − 2Δ/K) for K > K_c = 2Δ and R_ss = 0 below. Reference: Ott & Antonsen 2008, Chaos 18(3):037113.

Numerics¶

run(z0, n_steps) is the compute-kernel path: a tight RK4 loop on the real/imaginary components of z. This is dispatched across Rust / Mojo / Julia / Go / Python with bit-exact parity (scalar ODE, no reduction identities, no global sums — the only differences between backends are the rounding order of the k1..k4 accumulation, which matches exactly).

The scalar-output helpers K_c, steady_state_R and predict_from_oscillators stay native Python + optional Rust — they are O(1) arithmetic or O(N) percentile work and do not benefit from multi-language chains.

Classes¶

OAState `dataclass` ¶

OAState(z: complex, R: float, psi: float, K_c: float)

Ott-Antonsen mean-field state: order parameter and critical coupling.

OttAntonsenReduction ¶

OttAntonsenReduction(
    omega_0: float, delta: float, K: float, dt: float = 0.01
)

Ott-Antonsen mean-field reduction for globally-coupled Kuramoto.

The class stores (ω₀, Δ, K, dt) and exposes K_c, steady_state_R(), step(z), run(z0, n_steps) and predict_from_oscillators(omegas, K). run is dispatched across the 5-backend chain; the scalar helpers stay Python + optional Rust.

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

def __init__(
    self,
    omega_0: float,
    delta: float,
    K: float,
    dt: float = 0.01,
):
    omega_0 = _validate_finite_real(omega_0, name="omega_0")
    delta = _validate_finite_real(delta, name="delta")
    K = _validate_finite_real(K, name="K")
    dt = _validate_finite_real(dt, name="dt")
    if delta < 0:
        raise ValueError(f"delta (half-width) must be non-negative, got {delta}")
    if dt <= 0.0:
        raise ValueError(f"dt must be positive, got {dt}")
    self._omega_0 = omega_0
    self._delta = delta
    self._K = K
    self._dt = dt

Attributes¶

K_c `property` ¶

K_c: float

Critical coupling K_c = 2Δ.

Returns¶

float Critical coupling K_c = 2Δ.

Methods:¶

steady_state_R ¶

steady_state_R() -> float

Return the analytical steady-state R_ss = √(1 − 2Δ/K) for K > K_c.

Returns¶

float Return the analytical steady-state R_ss = √(1 − 2Δ/K) for K > K_c.

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

def steady_state_R(self) -> float:
    """Return the analytical steady-state ``R_ss = √(1 − 2Δ/K)`` for ``K > K_c``.

    Returns
    -------
    float
        Return the analytical steady-state ``R_ss = √(1 − 2Δ/K)`` for ``K > K_c``.
    """
    if _HAS_RUST_SCALAR:
        value = _rust_steady_state_r(self._delta, self._K)
    elif self.K_c >= self._K:
        value = 0.0
    else:
        value = (1.0 - 2.0 * self._delta / self._K) ** 0.5
    return _reduction_validation.validate_oa_steady_state_output(value)

step ¶

step(z: complex) -> complex

Single RK4 step on the OA ODE.

Parameters¶

z : complex Complex Ott-Antonsen order parameter.

Returns¶

complex The complex order parameter after one RK4 step.

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

def step(self, z: complex) -> complex:
    """Single RK4 step on the OA ODE.

    Parameters
    ----------
    z : complex
        Complex Ott-Antonsen order parameter.

    Returns
    -------
    complex
        The complex order parameter after one RK4 step.
    """
    z = _validate_finite_complex(z, name="z")
    re, im, _, _ = self._run_scalar(z.real, z.imag, n_steps=1)
    return complex(re, im)

run ¶

run(z0: complex, n_steps: int) -> OAState

Integrate n_steps RK4 steps; return the final OAState.

Parameters¶

z0 : complex Initial complex Ott-Antonsen order parameter. n_steps : int Number of integration steps to run.

Returns¶

OAState The final OAState after n_steps RK4 steps.

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

def run(self, z0: complex, n_steps: int) -> OAState:
    """Integrate ``n_steps`` RK4 steps; return the final ``OAState``.

    Parameters
    ----------
    z0 : complex
        Initial complex Ott-Antonsen order parameter.
    n_steps : int
        Number of integration steps to run.

    Returns
    -------
    OAState
        The final ``OAState`` after ``n_steps`` RK4 steps.
    """
    z0 = _validate_finite_complex(z0, name="z0")
    n_steps = _validate_positive_int(n_steps, name="n_steps")
    re, im, r, psi = self._run_scalar(z0.real, z0.imag, n_steps)
    return OAState(z=complex(re, im), R=r, psi=psi, K_c=self.K_c)

predict_from_oscillators ¶

predict_from_oscillators(
    omegas: FloatArray, K: float
) -> OAState

Fit a Lorentzian to omegas and return the relaxed OAState.

Parameters¶

omegas : FloatArray Natural frequencies in rad/s, shape (N,). K : float Global coupling strength.

Returns¶

OAState The relaxed OAState for the fitted Lorentzian.

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

def predict_from_oscillators(
    self,
    omegas: FloatArray,
    K: float,
) -> OAState:
    """Fit a Lorentzian to ``omegas`` and return the relaxed ``OAState``.

    Parameters
    ----------
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    K : float
        Global coupling strength.

    Returns
    -------
    OAState
        The relaxed ``OAState`` for the fitted Lorentzian.
    """
    omegas64 = _reduction_validation.validate_oa_frequency_sample(
        omegas,
        name="omegas",
    )
    K = _validate_finite_real(K, name="K")
    if _HAS_RUST_SCALAR:
        omega_0, delta = _rust_fit_lorentzian(omegas64)
    else:
        omega_0 = float(np.median(omegas64))
        q75, q25 = np.percentile(omegas64, [75, 25])
        delta = (q75 - q25) / 2.0 if q75 > q25 else 0.01
    reducer = OttAntonsenReduction(omega_0, delta, K, dt=self._dt)
    return reducer.run(complex(0.01, 0.0), n_steps=int(10.0 / self._dt))

Functions:¶

Dataclass: OAState¶

@dataclass
class OAState:
    z: complex   # Mean-field z = R·exp(iψ)
    R: float     # |z|, order parameter
    psi: float   # arg(z), mean phase
    K_c: float   # Critical coupling 2Δ

Rust Engine Functions¶

// RK4 integration of OA ODE, returns (z_re, z_im, R, psi)
pub fn oa_run(z_re, z_im, omega_0, delta, k_coupling, dt, n_steps) -> (f64, f64, f64, f64)

// Analytical steady-state R
pub fn steady_state_r_oa(delta, k_coupling) -> f64

// Fit Lorentzian (omega_0, delta) from frequency array
pub fn fit_lorentzian(omegas: &[f64]) -> (f64, f64)

Auto-Select Logic¶

All three functions (run, steady_state_R, predict_from_oscillators) use Rust when available. The RK4 loop in Rust is 38-96x faster because the Python path executes a Python-level loop with complex arithmetic per step.

7. Performance Benchmarks¶

Measured on Intel Core i5-11600K @ 3.90 GHz, 32 GB DDR4-2400. ω₀ = 0, Δ = 1, K = 4, z₀ = 0.01. Averaged over 100-1000 iterations.

Steps	Python (ms)	Rust (ms)	Speedup
100	0.203	0.0053	38.1x
500	1.836	0.0192	95.8x
2,000	7.183	0.0763	94.2x
10,000	20.811	0.3641	57.2x

Why 38-96x Speedup?¶

The OA reduction is a scalar complex ODE. Each RK4 step requires 4 evaluations of $f(z)$, each involving: - 1 complex multiply ($|z|^2$) - 2 complex multiply-adds - ~10 floating-point operations total

In Python, each step crosses the Python interpreter boundary for every arithmetic operation (complex + float overhead). In Rust, the entire loop runs in native code with no interpreter overhead. The 38-96x range reflects the interplay between Python interpreter overhead (dominant at high step counts → ~96x) and FFI call overhead (dominant at low step counts → ~38x). For the typical use case (500-2000 steps), the speedup is ~95x.

Memory Usage¶

Negligible: 2 floats (re, im) for state, 8 floats for RK4 stages. No heap allocation in the inner loop.

Test Coverage¶

Rust tests: 9 (reduction module in spo-engine)
Steady-state subcritical/supercritical, convergence, subcritical decay, rotation with ω₀, fit_lorentzian, zero steps, psi angle, critical coupling
Python tests: 12 (tests/test_ott_antonsen.py)
Critical coupling, below/above/at K_c, step finite, run convergence, run decay, negative delta raises, predict_from_oscillators, predict identical omegas, OAState fields, pipeline wiring
Source lines: 208 (Rust) + 99 (Python) = 307 total

8. Citations¶

Ott, E. & Antonsen, T. M. (2008). "Low dimensional behavior of large systems of globally coupled oscillators." Chaos 18(3):037113. DOI: 10.1063/1.2930766
Ott, E. & Antonsen, T. M. (2009). "Long time evolution of phase oscillator systems." Chaos 19(2):023117. DOI: 10.1063/1.3136851
Marvel, S. A., Mirollo, R. E., & Strogatz, S. H. (2009). "Identical phase oscillators with global sinusoidal coupling evolve by Möbius group action." Chaos 19(4):043104. DOI: 10.1063/1.3247089
Pikovsky, A. & Rosenblum, M. (2008). "Partially integrable dynamics of hierarchical populations of coupled oscillators." Physical Review Letters 101(26):264103. DOI: 10.1103/PhysRevLett.101.264103
Strogatz, S. H. (2000). "From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators." Physica D 143(1-4):1-20. DOI: 10.1016/S0167-2789(00)00094-4
Dörfler, F. & Bullo, F. (2014). "Synchronization in complex networks of phase oscillators: A survey." Automatica 50(6):1539-1564. DOI: 10.1016/j.automatica.2014.04.012
Kuramoto, Y. (1975). "Self-entrainment of a population of coupled non-linear oscillators." In Int. Symp. Mathematical Problems in Theoretical Physics, Lecture Notes in Physics 39:420-422. Springer.
Acebrón, J. A., Bonilla, L. L., Pérez Vicente, C. J., Ritort, F., & Spigler, R. (2005). "The Kuramoto model: A simple paradigm for synchronization phenomena." Reviews of Modern Physics 77:137-185. DOI: 10.1103/RevModPhys.77.137

Edge Cases and Limitations¶

|z₀| ≥ 1¶

The OA manifold requires $|z| \leq 1$. If seeded with $|z_0| > 1$, the cubic term $-|z|^2 z$ drives $|z|$ back below 1, but the transient may be non-physical. Always seed with $|z_0| \ll 1$ (e.g., 0.01).

Δ = 0 (Identical Oscillators)¶

When $\Delta = 0$, $K_c = 0$ — any positive coupling produces synchronisation. The steady-state $R_{ss} = \sqrt{1 - 0} = 1$. The OA ODE becomes $\dot{z} = -i\omega_0 z + \frac{K}{2}(z - |z|^2 z)$ — pure rotation plus cubic growth. Starting from $|z_0| < 1$, $|z| \to 1$ exponentially with rate $K/2$.

Non-Lorentzian Distributions¶

For Gaussian $g(\omega) = \frac{1}{\sigma\sqrt{2\pi}} e^{-(\omega-\omega_0)^2/2\sigma^2}$: - The OA ansatz is not exact - Using $\Delta = \sigma \cdot \sqrt{2\ln 2}$ (FWHM matching) gives $R_{ss}$ accurate to $\sim 5\%$ for moderate coupling - For bimodal distributions, the OA ansatz does not apply

Numerical Stability¶

The RK4 integrator is unconditionally stable for $K \Delta t < 4$. For typical parameters ($K \leq 10$, $\Delta t = 0.01$), this is always satisfied. If $K > 100$, reduce $\Delta t$ accordingly.

Finite-$N$ Corrections¶

For $N$ oscillators (not $N \to \infty$): - $R$ fluctuates around $R_{ss}$ with variance $\sim 1/N$ - The transition at $K_c$ is rounded: $R > 0$ for all $K > 0$ (finite-size precursors) - The predict_from_oscillators method accounts for the Lorentzian fitting uncertainty but not the finite-$N$ fluctuations

Derivation of the Ott-Antonsen Ansatz¶

Step 1: Continuity Equation¶

For $N \to \infty$ oscillators with density $f(\theta, \omega, t)$ on $S^1$, the continuity equation reads:

\[\frac{\partial f}{\partial t} + \frac{\partial}{\partial \theta}\big(f \cdot v(\theta, \omega, t)\big) = 0\]

where $v = \omega + K R \sin(\psi - \theta)$ is the velocity field from the Kuramoto model with mean-field coupling.

Step 2: Fourier Expansion¶

The density admits a Fourier series in $\theta$:

\[f(\theta, \omega, t) = \frac{1}{2\pi}\left[1 + \sum_{n=1}^{\infty} \hat{f}_n(\omega, t) e^{in\theta} + \text{c.c.}\right]\]

Substituting into the continuity equation yields an infinite hierarchy of ODEs for the Fourier coefficients $\hat{f}_n$.

Step 3: The Ansatz¶

Ott and Antonsen (2008) observed that if we restrict to the manifold:

\[\hat{f}_n(\omega, t) = \alpha(\omega, t)^n\]

then the infinite hierarchy collapses to a single ODE for $\alpha$:

\[\dot{\alpha} = -i\omega\alpha + \frac{K}{2}(\bar{z}\alpha^2 - z)\]

where $z = R e^{i\psi}$ is the Kuramoto order parameter, defined by:

\[z(t) = \int_{-\infty}^{\infty} \int_0^{2\pi} e^{i\theta} f(\theta, \omega, t)\,d\theta\,d\omega\]

Step 4: Lorentzian Closure¶

For $g(\omega) = \frac{\Delta/\pi}{(\omega - \omega_0)^2 + \Delta^2}$, the integral over $\omega$ can be evaluated by closing the contour in the lower half-plane, picking up the single pole at $\omega = \omega_0 - i\Delta$:

\[z(t) = \overline{\alpha(\omega_0 - i\Delta, t)}\]

This yields the closed ODE:

\[\dot{z} = -(\Delta + i\omega_0)z + \frac{K}{2}(z - |z|^2 z)\]

The key insight is that the Lorentzian distribution has a single pole in the lower half-plane, which makes the contour integral exact. Distributions with multiple poles (e.g., bimodal) or no poles (e.g., Gaussian, which has an essential singularity) break the closure.

Rust Implementation Details¶

oa_deriv — Inner Loop Kernel¶

The derivative function oa_deriv computes both real and imaginary parts of $f(z)$ without heap allocation:

fn oa_deriv(re: f64, im: f64, omega_0: f64, delta: f64, half_k: f64) -> (f64, f64) {
    let abs_sq = re * re + im * im;
    let lin_re = -delta * re + omega_0 * im;
    let lin_im = -delta * im - omega_0 * re;
    let cubic_factor = half_k * (1.0 - abs_sq);
    (lin_re + cubic_factor * re, lin_im + cubic_factor * im)
}

Each RK4 step calls oa_deriv 4 times (k1-k4), for a total of ~40 floating-point operations per step. The Rust compiler inlines this function and uses SSE2/AVX for the multiply-adds.

fit_lorentzian — Robust Estimation¶

The Lorentzian fitting uses order statistics rather than maximum-likelihood estimation: - Median for $\omega_0$: breakdown point 50%, robust to outliers - IQR/2 for $\Delta$: exploits the exact relationship $\text{IQR} = 2\Delta$ for a Cauchy distribution

The sort is $O(n \log n)$, but $n$ is the number of oscillators (typically $\leq 10{,}000$), so the cost is negligible compared to the integration.

FFI Transfer¶

The oa_run function returns a tuple (f64, f64, f64, f64), which PyO3 converts to a Python tuple with zero-copy semantics (scalar values only, no array allocation). The fit_lorentzian function accepts a PyReadonlyArray1<f64> reference, avoiding any data copy from NumPy to Rust.

Troubleshooting¶

Issue: R Does Not Converge¶

Symptom: After many steps, $R$ oscillates or drifts instead of settling to $R_{ss}$.

Diagnosis: Check whether $\omega_0 \neq 0$. In the rotating frame, $R$ converges while $\psi$ rotates at frequency $\omega_0$. The complex $z$ traces a circle of radius $R_{ss}$, so $R = |z|$ converges even though $z$ itself does not settle to a fixed point.

If $R$ truly oscillates, verify that $K \Delta t < 4$ (RK4 stability bound). For $K = 100$ and $\Delta t = 0.01$, $K \Delta t = 1.0$ is safe.

Issue: predict_from_oscillators Returns R ≈ 0 for Clearly Coupled System¶

Diagnosis: The Lorentzian fit via IQR assumes unimodal, heavy-tailed $g(\omega)$. For bimodal or Gaussian distributions, the estimated $\Delta$ may be too large, yielding $K_c > K$.

Solution: Manually set $(\omega_0, \Delta)$ from known distribution parameters rather than using the automatic fitting.

Issue: R_ss = 0 Despite Supercritical K¶

Diagnosis: The steady_state_R function requires K > 2Δ. If you pass K as the total coupling but the per-pair coupling is $K/N$, you need to use the per-pair value multiplied by $N$ (i.e., the effective global coupling strength).

Issue: Mismatch Between OA and Full Simulation¶

Diagnosis: Common causes: 1. Non-Lorentzian frequency distribution (OA is approximate) 2. Finite-$N$ fluctuations ($R$ deviates by $\sim 1/\sqrt{N}$) 3. Phase frustration ($\alpha_{ij} \neq 0$) — OA assumes zero frustration 4. Non-global coupling (sparse $K_{ij}$) — OA assumes all-to-all

For cases 1-2, the mismatch should decrease with larger $N$ and Lorentzian-like distributions. For cases 3-4, the OA reduction is not applicable.

Integration with Other SPO Modules¶

With UPDEEngine¶

The OttAntonsenReduction validates full simulations:

# Quick check: will this parameter set synchronise?
oa = OttAntonsenReduction(omega_0=0, delta=0.5, K=3.0)
if oa.steady_state_R() > 0.5:
    # Worth running the expensive N-oscillator simulation
    engine = UPDEEngine(n_oscillators=256, dt=0.01)
    # ... run full sim ...

With SSGF Costs¶

The predicted $R_{ss}$ can initialise the SSGF sync deficit cost:

\[C_1 = 1 - R \approx 1 - R_{ss}\]

This provides a baseline expectation for the geometry optimisation. If the full simulation yields $R < R_{ss}$, the geometry is suboptimal and the SSGF engine should increase the learning rate.

With RegimeManager¶

The OA reduction enables fast regime classification:

$K / K_c$	Regime	$R_{ss}$
$< 0.5$	Incoherent	0
$0.5 - 1.0$	Subcritical	0
$1.0 - 1.5$	Onset	$0 - 0.58$
$1.5 - 3.0$	Partial sync	$0.58 - 0.82$
$> 3.0$	Full sync	$> 0.82$

The supervisor can use these thresholds to decide whether to engage geometry control or whether the current coupling is sufficient.

\(K / K_c\)	Regime	\(R_{ss}\)
\(< 0.5\)	Incoherent	0
\(0.5 - 1.0\)	Subcritical	0
\(1.0 - 1.5\)	Onset	\(0 - 0.58\)
\(1.5 - 3.0\)	Partial sync	\(0.58 - 0.82\)
\(> 3.0\)	Full sync	\(> 0.82\)

Field	Type	Range	Meaning
`z`	`complex`	$	z
`R`	`float`	\([0, 1]\)	$
`psi`	`float`	\((-\pi, \pi]\)	\(\arg(z)\), mean phase
`K_c`	`float`	\(\geq 0\)	Critical coupling \(2\Delta\)

Ott-Antonsen Mean-Field Reduction¶

1. Mathematical Formalism¶

The Ott-Antonsen Ansatz¶

Decomposition of the ODE¶

Steady-State Solution¶

Critical Coupling¶

Integration Method¶

2. Theoretical Context¶

Why Mean-Field Reduction?¶

Domain of Validity¶

Historical Context¶

Applications in SPO¶

Lorentzian Fitting¶

3. Pipeline Position¶

Input Contracts¶

Output Contract¶

Accelerator Output Contract¶

4. Features¶

5. Usage Examples¶

Basic: Predict R for Given Parameters¶

Fast Diagnostic from Oscillator Frequencies¶

Parameter Sweep¶

Validation Against Full Simulation¶

6. Technical Reference¶

Class: OttAntonsenReduction¶

reduction ¶

Dynamics¶

Numerics¶

Classes¶

OAState dataclass ¶

OttAntonsenReduction ¶

Attributes¶

K_c property ¶

Returns¶

Methods:¶

steady_state_R ¶

Returns¶

step ¶

Parameters¶

Returns¶

run ¶

Parameters¶

Returns¶

predict_from_oscillators ¶

Parameters¶

Returns¶

Functions:¶

Dataclass: OAState¶

Rust Engine Functions¶

Auto-Select Logic¶

7. Performance Benchmarks¶

Why 38-96x Speedup?¶

Memory Usage¶

Test Coverage¶

8. Citations¶

Edge Cases and Limitations¶

|z₀| ≥ 1¶

Δ = 0 (Identical Oscillators)¶

Non-Lorentzian Distributions¶

Numerical Stability¶

Finite-\(N\) Corrections¶

Derivation of the Ott-Antonsen Ansatz¶

Step 1: Continuity Equation¶

Step 2: Fourier Expansion¶

Step 3: The Ansatz¶

Step 4: Lorentzian Closure¶

Rust Implementation Details¶

oa_deriv — Inner Loop Kernel¶

fit_lorentzian — Robust Estimation¶

FFI Transfer¶

Troubleshooting¶

Issue: R Does Not Converge¶

Issue: predict_from_oscillators Returns R ≈ 0 for Clearly Coupled System¶

Issue: R_ss = 0 Despite Supercritical K¶

Issue: Mismatch Between OA and Full Simulation¶

Integration with Other SPO Modules¶

With UPDEEngine¶

With SSGF Costs¶

With RegimeManager¶

OAState `dataclass` ¶

K_c `property` ¶