Hypergraph (k-Body) Kuramoto Engine

1. Mathematical Formalism

Generalised k-Body Coupling

The hypergraph Kuramoto model extends pairwise coupling to arbitrary higher-order interactions. For a \(k\)-hyperedge \(e = \{i_1, \ldots, i_k\}\) with coupling strength \(\sigma_e\), the contribution to oscillator \(i_m\) in the hyperedge is:

\[\dot{\theta}_{i_m}\big|_e = \sigma_e \cdot \sin\!\left(\sum_{\substack{j \in e \\ j \neq m}} \theta_{i_j} - (k-1)\,\theta_{i_m}\right)\]

Using the total phase sum \(\Phi_e = \sum_{j \in e} \theta_{i_j}\), the argument simplifies to:

\[\dot{\theta}_{i_m}\big|_e = \sigma_e \cdot \sin(\Phi_e - k\,\theta_{i_m})\]

Full Dynamics

The complete equation of motion for oscillator \(i\) combines natural frequency, optional pairwise coupling, all hyperedge contributions, and external drive:

\[\frac{d\theta_i}{dt} = \omega_i + \underbrace{\sum_{j=1}^{N} K_{ij} \sin(\theta_j - \theta_i - \alpha_{ij})}_{\text{pairwise (optional)}} + \underbrace{\sum_{e \ni i} \sigma_e \sin(\Phi_e - |e|\,\theta_i)}_{\text{hypergraph k-body}} + \underbrace{\zeta \sin(\Psi - \theta_i)}_{\text{external drive}}\]

Special Cases

Order \(k\)	Argument	Model
\(k = 2\)	\(\sin(\theta_j - \theta_i)\)	Standard Kuramoto
\(k = 3\)	\(\sin(\theta_j + \theta_k - 2\theta_i)\)	Simplicial (Tanaka & Aoyagi 2011)
\(k = 4\)	\(\sin(\theta_j + \theta_k + \theta_l - 3\theta_i)\)	4-body (Bick et al. 2023)
Mixed	Combination of above	This engine

The HypergraphEngine supports mixed-order interactions: some edges can be pairwise, some 3-body, some 4-body, in the same system.

Integration Method

Euler forward integration:

\[\theta_i(t + \Delta t) = \left(\theta_i(t) + \Delta t \cdot \dot{\theta}_i(t)\right) \bmod 2\pi\]

The run() method executes the full loop in Rust via spo_kernel when available, with flat-encoded hyperedges for zero-copy FFI transfer.

---

2. Theoretical Context

Why Hypergraph Coupling?

Real complex systems exhibit interactions that cannot be reduced to pairwise connections:

• Neural ensembles: Synaptic transmission involves presynaptic calcium dynamics, vesicle release probability, and postsynaptic receptor occupancy — a fundamentally multi-party process (Petri et al. 2014).
• Gene regulatory networks: Transcription factor complexes require \(k\) proteins to co-localise, creating \(k\)-body logical gates (Battiston et al. 2020).
• Social contagion: Adoption of behaviours requires reinforcement from multiple peers simultaneously (threshold models), not just pairwise influence (Iacopini et al. 2019).
• Ecological webs: Pollinator-plant-herbivore triangles exhibit dynamics absent from pairwise predator-prey models.

Theoretical Predictions

Tanaka & Aoyagi (2011) showed that \(k\)-body coupling with \(k \geq 3\) produces multistable attractors — multiple distinct synchronised states coexist at the same coupling strength. This is qualitatively different from pairwise Kuramoto, which has at most one stable synchronised state for given parameters.

Skardal & Arenas (2019) proved that higher-order coupling induces explosive (first-order) synchronisation transitions, where the order parameter \(R\) jumps discontinuously. The critical coupling scales as:

\[K_c^{(k)} \propto \frac{\Delta}{\langle R^{k-2} \rangle}\]

where \(\Delta\) is the frequency spread and \(\langle R^{k-2} \rangle\) is the \((k-2)\)-th moment of the order parameter distribution.

Topology and Algebraic Structure

A hypergraph \(H = (V, E)\) consists of a vertex set \(V\) and a collection of hyperedges \(E\) where each \(e \in E\) is a subset \(e \subseteq V\) with \(|e| \geq 2\). Unlike simplicial complexes, hypergraphs do not require the downward closure property: a 3-body edge \(\{a, b, c\}\) can exist without the pairwise edges \(\{a,b\}\), \(\{a,c\}\), \(\{b,c\}\).

The adjacency tensor generalises the adjacency matrix:

\[A^{(k)}_{i_1 \ldots i_k} = \begin{cases} \sigma_e & \text{if } \{i_1, \ldots, i_k\} = e \in E \\ 0 & \text{otherwise} \end{cases}\]

For \(k = 2\), this reduces to the standard adjacency matrix \(A_{ij}\).

Ott-Antonsen Reduction for Hypergraphs

For all-to-all \(k\)-body coupling with Lorentzian frequency distribution, the Ott-Antonsen ansatz yields a reduced equation for the mean-field \(z = R e^{i\psi}\) (Skardal & Arenas 2020):

\[\dot{z} = -(\Delta + i\omega_0) z + \frac{\sigma_k}{2} \left(z^{k-1} \bar{z}^{k-2} - |z|^{2(k-1)} z\right)\]

This generalises the standard OA equation (\(k = 2\)) and explains why higher \(k\) produces sharper transitions: the effective nonlinearity increases with \(k\).

Applications in Neuroscience

Recent work by Petri et al. (2014) and Reimann et al. (2017) showed that brain functional networks exhibit significant higher-order structure (simplicial complexes of dimension up to 7 in the Blue Brain connectome). The hypergraph engine enables simulation of these structures with explicit control over which groups of neurons participate in collective interactions.

Key parameters for neural hypergraphs: - Order distribution: Which \(k\) values are present (typically \(k = 2, 3, 4\)) - Strength hierarchy: \(\sigma_3 < \sigma_2\) (higher-order weaker) - Spatial constraint: Hyperedges limited to local circuit motifs

Comparison with SimplicialEngine

SimplicialEngine is a specialised, optimised version for the \(k=3\) case. It uses the trig factorisation \(2 S_i C_i\) to avoid looping over edges. HypergraphEngine is general: it handles any mix of orders but loops over each hyperedge explicitly. Use SimplicialEngine when all higher-order interactions are exactly 3-body for better performance.

Feature	SimplicialEngine	HypergraphEngine
Orders supported	2 + 3	Any \(k \geq 2\)
Edge representation	Implicit (all-to-all 3-body)	Explicit edge list
3-body complexity	\(O(N^2)\) via factorisation	$O(
Mixed orders	No	Yes
Pairwise via K_nm	Yes	Optional

---

3. Pipeline Position

Oscillators.extract() ──→ θ, ω
                               │
CouplingBuilder.build() ──→ K_nm, α (optional pairwise)
                               │
Hyperedge definitions ──→ [(nodes, σ), ...]
                               │
                               ↓
     ┌──── HypergraphEngine(n, dt, hyperedges) ────┐
     │                                              │
     │  Input:  θ, ω, [optional K_nm, α, ζ, Ψ]    │
     │  Edges:  list of Hyperedge(nodes, strength)  │
     │  Method: Euler (Rust FFI or Python)          │
     │  Output: θ_new ∈ [0, 2π)^N                  │
     │                                              │
     └──────────────────────────────────────────────┘
                               │
                               ↓
              compute_order_parameter(θ_new) → R, ψ

FFI Encoding

Hyperedges are serialised for Rust FFI as three flat arrays: - edge_nodes: int64[] — concatenated node indices for all edges - edge_offsets: int64[] — start index in edge_nodes for each edge - edge_strengths: float64[] — coupling strength per edge

Direct Go, Julia, and Mojo accelerator wrappers validate this encoding before loading optional runtimes: offsets must start at zero and increase strictly, each encoded hyperedge must contain at least two distinct valid oscillator indices, strengths and state vectors must be finite real float64 buffers, and optional pairwise K_nm/alpha flat buffers must be empty or exactly N*N with zero diagonals. Zero-step direct calls return the torus-normalised input phases without requiring the optional runtime.

The public HypergraphEngine.run() dispatcher and Rust wrapper replay the same output validator before publishing optional-backend results. Returned phase vectors must keep oscillator cardinality, contain only finite real values, and remain on the torus in [0, 2*pi); loader/runtime unavailability can still fall back to Python, but malformed backend physics evidence raises immediately. Public phase/frequency arrays, optional pairwise and phase-lag matrices, hyperedge strengths, scalar controls, step counts, direct backend vectors and index buffers, optional-backend outputs, and direct Julia raw returns reject numeric-string aliases before float or integer coercion.

This avoids Python object overhead and enables zero-copy transfer.

Input Contracts

Parameter	Type	Shape	Source
phases	NDArray[float64]	(N,)	Previous step
omegas	NDArray[float64]	(N,)	Oscillators
pairwise_knm	NDArray[float64] \| None	(N, N)	CouplingBuilder (optional)
alpha	NDArray[float64] \| None	(N, N)	Phase-lag (optional)
zeta	float	scalar	External drive
psi	float	scalar	Drive phase

---

4. Features

• Arbitrary \(k\)-body coupling — any order \(k \geq 2\)
• Mixed-order edges — combine pairwise, 3-body, 4-body in one system
• Dynamic edge addition — add_edge() and add_all_to_all(order)
• Optional pairwise coupling via standard \(K_{nm}\) matrix
• External drive (\(\zeta\), \(\Psi\)) for entrainment
• Full Rust FFI acceleration with flat-encoded edge transfer
• Direct Go/Julia/Mojo boundary validation before optional runtime loading
• Public optional-backend output validation before phase vectors publish
• Pre-coercion numeric-string rejection on public and direct boundaries
• Order parameter computation — built-in order_parameter() method
• step() and run() — single-step or batch integration

---

5. Usage Examples

Basic: 3-Body All-to-All

import numpy as np
from scpn_phase_orchestrator.upde.hypergraph import HypergraphEngine, Hyperedge
from scpn_phase_orchestrator.upde.order_params import compute_order_parameter

N = 16
engine = HypergraphEngine(N, dt=0.01)
engine.add_all_to_all(order=3, strength=0.5)

print(f"Number of 3-body edges: {engine.n_edges}")  # C(16,3) = 560

phases = np.random.default_rng(42).uniform(0, 2 * np.pi, N)
omegas = np.ones(N)

phases_final = engine.run(phases, omegas, n_steps=1000)
R, _ = compute_order_parameter(phases_final)
print(f"R = {R:.4f}")

Mixed-Order: Pairwise + 3-Body + 4-Body

import numpy as np
from scpn_phase_orchestrator.upde.hypergraph import HypergraphEngine

N = 8
engine = HypergraphEngine(N, dt=0.01)

# Pairwise edges (specific pairs, not K_nm matrix)
engine.add_edge((0, 1), strength=1.0)
engine.add_edge((2, 3), strength=1.0)
engine.add_edge((4, 5), strength=1.0)

# 3-body triangles
engine.add_edge((0, 1, 2), strength=0.5)
engine.add_edge((3, 4, 5), strength=0.5)

# 4-body interaction linking both groups
engine.add_edge((1, 2, 4, 5), strength=0.3)

phases = np.zeros(N)
phases[4:] = np.pi  # Two anti-phase clusters
omegas = np.ones(N)

phases_final = engine.run(phases, omegas, n_steps=2000)
R = engine.order_parameter(phases_final)
print(f"R = {R:.4f}")

Combined with Pairwise K_nm Matrix

import numpy as np
from scpn_phase_orchestrator.upde.hypergraph import HypergraphEngine
from scpn_phase_orchestrator.coupling.knm import CouplingBuilder

N = 12
builder = CouplingBuilder(N, method="exponential", K_base=0.5, decay=0.3)
knm, alpha = builder.build()

engine = HypergraphEngine(N, dt=0.01)
engine.add_all_to_all(order=3, strength=0.2)

phases = np.random.default_rng(42).uniform(0, 2 * np.pi, N)
omegas = np.linspace(0.8, 1.2, N)

# Pairwise via K_nm + 3-body via hyperedges
phases_final = engine.run(
    phases, omegas, n_steps=1000,
    pairwise_knm=knm, alpha=alpha,
)

Measuring Phase Transitions

import numpy as np
from scpn_phase_orchestrator.upde.hypergraph import HypergraphEngine

N = 20
omegas = np.random.default_rng(42).standard_cauchy(N) * 0.05
phases_init = np.random.default_rng(42).uniform(0, 2 * np.pi, N)

sigma_values = np.linspace(0, 2, 40)
R_values = []

for sigma in sigma_values:
    eng = HypergraphEngine(N, dt=0.01)
    eng.add_all_to_all(order=3, strength=sigma)
    p = eng.run(phases_init.copy(), omegas, n_steps=2000)
    R_values.append(eng.order_parameter(p))

# Find critical coupling
for i, (s, r) in enumerate(zip(sigma_values, R_values)):
    if r > 0.5:
        print(f"σ_c ≈ {s:.2f} (R = {r:.3f})")
        break

---

6. Technical Reference

Class: HypergraphEngine

hypergraph

Hypergraph Kuramoto with arbitrary k-body interactions beyond pairwise.

Exposes a 5-backend fallback chain.

Extends the standard Kuramoto model with k-body coupling terms for any k ≥ 2. The standard model (k=2) and simplicial model (k=3) are special cases.

For a k-hyperedge {i₁, …, iₖ}, the coupling on oscillator iₘ is

σₖ · sin( Σ_{j≠m} θ_{iⱼ} − (k−1)·θ_{iₘ} )

which generalises

k = 2:  sin(θ_j − θ_i)            — standard Kuramoto
k = 3:  sin(θ_j + θ_k − 2·θ_i)    — simplicial / triadic

The engine also accepts a dense pairwise coupling matrix pairwise_knm (optional) and an external-drive field (ζ, ψ).

Numerics

The pairwise-derivative loop uses the Rust kernel's sin(θ_j − θ_i) = sin(θ_j)·cos(θ_i) − cos(θ_j)·sin(θ_i) expansion in the alpha == 0 fast path so that floating-point rounding matches Rust (spo-engine/src/hypergraph.rs) bit-for-bit. Alpha ≠ 0 falls back to the direct sin(θ_j − θ_i − α) form in all five backends.

References

Tanaka & Aoyagi 2011, Phys. Rev. Lett. 106:224101. Skardal & Arenas 2019, Comm. Phys. 2:22. Bick, Gross, Harrington & Schaub 2023, Nat. Rev. Physics 5:307-317.

Classes

Hyperedge dataclass

Hyperedge(nodes: tuple[int, ...], strength: float = 1.0)

A k-body interaction among oscillators.

Attributes

nodes: Tuple of oscillator indices in this hyperedge.
strength: Coupling strength σₖ for this hyperedge.

Attributes

order property

order: int

Return the number of oscillators participating in the hyperedge.

Returns

int Return the number of oscillators participating in the hyperedge.

HypergraphEngine

HypergraphEngine(
    n_oscillators: int,
    dt: float,
    hyperedges: list[Hyperedge] | None = None,
)

Kuramoto engine with arbitrary k-body hypergraph coupling.

Supports mixed-order interactions: some edges can be pairwise, some 3-body, some 4-body, etc. Each Hyperedge specifies which oscillators participate and the coupling strength.

Initialise a validated hypergraph Kuramoto engine.

Parameters

n_oscillators : int Positive number of oscillators in the simulated system. dt : float Positive finite explicit-Euler step size. hyperedges : list[Hyperedge] | None Optional initial hyperedge definitions to validate and store.

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

def __init__(
    self,
    n_oscillators: int,
    dt: float,
    hyperedges: list[Hyperedge] | None = None,
) -> None:
    """Initialise a validated hypergraph Kuramoto engine.

    Parameters
    ----------
    n_oscillators : int
        Positive number of oscillators in the simulated system.
    dt : float
        Positive finite explicit-Euler step size.
    hyperedges : list[Hyperedge] | None
        Optional initial hyperedge definitions to validate and store.
    """
    self._n = _validate_positive_int(n_oscillators, name="n_oscillators")
    self._dt = _validate_positive_float(dt, name="dt")
    self._hyperedges = [
        _validate_hyperedge(edge, n_oscillators=self._n)
        for edge in (hyperedges or [])
    ]

Attributes

n_edges property

n_edges: int

Return the number of configured hyperedges.

Returns

int Return the number of configured hyperedges.

Methods:

add_edge

add_edge(
    nodes: tuple[int, ...], strength: float = 1.0
) -> None

Validate and append one explicit k-body hyperedge.

Parameters

nodes : tuple[int, ...] Indices of the oscillators participating in the hyperedge. strength : float Coupling strength assigned to the hyperedge(s).

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

def add_edge(self, nodes: tuple[int, ...], strength: float = 1.0) -> None:
    """Validate and append one explicit k-body hyperedge.

    Parameters
    ----------
    nodes : tuple[int, ...]
        Indices of the oscillators participating in the hyperedge.
    strength : float
        Coupling strength assigned to the hyperedge(s).
    """
    edge = _validate_hyperedge(
        Hyperedge(nodes=nodes, strength=strength),
        n_oscillators=self._n,
    )
    self._hyperedges.append(edge)

add_all_to_all

add_all_to_all(order: int, strength: float = 1.0) -> None

Add all C(N, order) hyperedges of given order.

Parameters

order : int Interaction order (number of oscillators per hyperedge). strength : float Coupling strength assigned to the hyperedge(s).

Raises

ValueError If order is outside 2..N.

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

def add_all_to_all(self, order: int, strength: float = 1.0) -> None:
    """Add all C(N, order) hyperedges of given order.

    Parameters
    ----------
    order : int
        Interaction order (number of oscillators per hyperedge).
    strength : float
        Coupling strength assigned to the hyperedge(s).

    Raises
    ------
    ValueError
        If ``order`` is outside ``2..N``.
    """
    from itertools import combinations

    order = _validate_positive_int(order, name="order")
    if order > self._n:
        raise ValueError(f"order must be <= n_oscillators, got {order!r}")
    validated_edge = _validate_hyperedge(
        Hyperedge(nodes=tuple(range(order)), strength=strength),
        n_oscillators=self._n,
    )
    for combo in combinations(range(self._n), order):
        self._hyperedges.append(
            Hyperedge(nodes=combo, strength=validated_edge.strength)
        )

step

step(
    phases: FloatArray,
    omegas: FloatArray,
    pairwise_knm: FloatArray | None = None,
    alpha: FloatArray | None = None,
    zeta: float = 0.0,
    psi: float = 0.0,
) -> FloatArray

One explicit-Euler step.

Parameters

phases : FloatArray Oscillator phases in radians, shape (N,). omegas : FloatArray Natural frequencies in rad/s, shape (N,). pairwise_knm : FloatArray | None Optional pairwise coupling matrix (N, N), or None. alpha : FloatArray | None Phase-lag matrix in radians, shape (N, N), or None for no lag. zeta : float External drive strength ζ. psi : float External drive reference phase Ψ in radians.

Returns

FloatArray The phases after one explicit-Euler hypergraph step.

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

def step(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    pairwise_knm: FloatArray | None = None,
    alpha: FloatArray | None = None,
    zeta: float = 0.0,
    psi: float = 0.0,
) -> FloatArray:
    """One explicit-Euler step.

    Parameters
    ----------
    phases : FloatArray
        Oscillator phases in radians, shape ``(N,)``.
    omegas : FloatArray
        Natural frequencies in rad/s, shape ``(N,)``.
    pairwise_knm : FloatArray | None
        Optional pairwise coupling matrix ``(N, N)``, or ``None``.
    alpha : FloatArray | None
        Phase-lag matrix in radians, shape ``(N, N)``, or ``None`` for no lag.
    zeta : float
        External drive strength ``ζ``.
    psi : float
        External drive reference phase ``Ψ`` in radians.

    Returns
    -------
    FloatArray
        The phases after one explicit-Euler hypergraph step.
    """
    return self.run(
        phases,
        omegas,
        n_steps=1,
        pairwise_knm=pairwise_knm,
        alpha=alpha,
        zeta=zeta,
        psi=psi,
    )

run

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

Integrate n_steps Euler steps via the fastest backend; return phases.

Parameters

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

Returns

FloatArray The final phases after n_steps hypergraph steps.

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

def run(
    self,
    phases: FloatArray,
    omegas: FloatArray,
    n_steps: int,
    pairwise_knm: FloatArray | None = None,
    alpha: FloatArray | None = None,
    zeta: float = 0.0,
    psi: float = 0.0,
) -> FloatArray:
    """Integrate ``n_steps`` Euler steps via the fastest backend; return phases.

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

    Returns
    -------
    FloatArray
        The final phases after ``n_steps`` hypergraph 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,))
    knm_flat = _validate_optional_state_array(
        pairwise_knm,
        name="pairwise_knm",
        shape=(self._n, self._n),
    )
    alpha_flat = _validate_optional_state_array(
        alpha,
        name="alpha",
        shape=(self._n, self._n),
    )
    zeta = _validate_finite_float(zeta, name="zeta")
    psi = _validate_finite_float(psi, name="psi")
    en, eo, es = self._encode_edges()
    backend_fn = _dispatch()
    if backend_fn is not None:
        return _validate_backend_output(
            backend_fn(
                phases64,
                omegas64,
                self._n,
                en,
                eo,
                es,
                knm_flat,
                alpha_flat,
                float(zeta),
                float(psi),
                float(self._dt),
                int(n_steps),
            ),
            n=self._n,
        )
    return _validate_backend_output(
        _python_run(
            phases64,
            omegas64,
            self._n,
            en,
            eo,
            es,
            knm_flat,
            alpha_flat,
            float(zeta),
            float(psi),
            float(self._dt),
            int(n_steps),
        ),
        n=self._n,
    )

order_parameter

order_parameter(phases: FloatArray) -> float

Compute the standard Kuramoto R = ||.

Parameters

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

Returns

float The Kuramoto order parameter R.

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

def order_parameter(self, phases: FloatArray) -> float:
    """Compute the standard Kuramoto R = |<exp(iθ)>|.

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

    Returns
    -------
    float
        The Kuramoto order parameter ``R``.
    """
    phases64 = _validate_state_array(
        phases,
        name="phases",
        shape=(self._n,),
    )
    return float(np.abs(np.mean(np.exp(1j * phases64))))

Functions:

Dataclass: Hyperedge

@dataclass
class Hyperedge:
    nodes: tuple[int, ...]   # oscillator indices in this hyperedge
    strength: float = 1.0    # coupling strength σ_e

    @property
    def order(self) -> int:  # k = number of nodes
        return len(self.nodes)

Rust Engine Function

pub fn hypergraph_run(
    phases: &[f64],          // length N
    omegas: &[f64],          // length N
    n: usize,
    edges: &[Hyperedge],     // reconstructed from flat encoding
    pairwise_knm: &[f64],   // length N*N or empty
    alpha: &[f64],           // length N*N or empty
    zeta: f64,
    psi: f64,
    dt: f64,
    n_steps: usize,
) -> Vec<f64>                // length N

FFI Flat Encoding

The Python wrapper serialises hyperedges as:

edge_nodes = np.array([0, 1, 2, 1, 2, 3], dtype=np.int64)  # two 3-body edges
edge_offsets = np.array([0, 3], dtype=np.int64)               # starts at 0 and 3
edge_strengths = np.array([0.5, 0.5], dtype=np.float64)

Rust reconstructs Hyperedge structs from these arrays with zero allocation overhead per edge beyond the initial Vec::with_capacity.

Auto-Select Logic

Rust path is used only when _HAS_RUST is True and there are hyperedges (self._hyperedges is non-empty). If only pairwise coupling is used (no hyperedges), the Python loop is trivially fast and the Rust overhead is not warranted.

---

7. Performance Benchmarks

Measured on Intel Core i5-11600K @ 3.90 GHz, 32 GB DDR4-2400. 200 Euler steps, all-to-all 3-body edges (\(C(N,3)\) hyperedges), \(\sigma = 0.5\). Averaged over 3-20 iterations.

N	Edges \(C(N,3)\)	Python (ms)	Rust (ms)	Speedup
8	56	52.85	0.52	101.7x
16	560	652.20	7.92	82.4x
24	2,024	2,152.21	30.75	70.0x

Scaling Analysis

The computation is \(O(|E| \cdot k \cdot T)\) where \(|E|\) is the number of hyperedges, \(k\) is the average order, and \(T\) is the number of steps. For all-to-all 3-body edges, \(|E| = C(N,3) = O(N^3/6)\), making the total cost \(O(N^3 \cdot T)\).

The Rust speedup is dramatic (70-100x) because: 1. Python per-edge loop overhead is eliminated 2. No Python object allocation per edge per step 3. Cache-friendly flat array access in Rust

For large \(N\), consider using SimplicialEngine instead (all-to-all 3-body via trig factorisation at \(O(N^2)\), not \(O(N^3)\)).

Memory Usage

Rust: \(O(N + |E| \cdot k)\) for edge storage + \(O(N)\) working arrays. Python: Same edge storage + NumPy temporaries for pairwise terms.

Test Coverage

• Rust tests: 9 (hypergraph module in spo-engine)
• Free rotation, pairwise-only synchronisation, 3-body effect, 4-body, mixed-order, external drive, zero steps, order parameter synchronised, order parameter uniform
• Python tests: 11 (tests/test_hypergraph.py)
• Source lines: 242 (Rust) + 176 (Python) = 418 total

---

8. Citations

1. Tanaka, T. & Aoyagi, T. (2011). "Multistable attractors in a network of phase oscillators with three-body interactions." Physical Review Letters 106:224101. DOI: 10.1103/PhysRevLett.106.224101

2. Skardal, P. S. & Arenas, A. (2019). "Abrupt desynchronization and extensive multistability in globally coupled oscillator simplexes." Physical Review Letters 122:248301. DOI: 10.1103/PhysRevLett.122.248301

3. Bick, C., Gross, E., Harrington, H. A., & Schaub, M. T. (2023). "What are higher-order networks?" Nature Reviews Physics 5:307-317. DOI: 10.1038/s42254-023-00573-y

4. Battiston, F., Cencetti, G., Iacopini, I., et al. (2020). "Networks beyond pairwise interactions: Structure and dynamics." Physics Reports 874:1-92. DOI: 10.1016/j.physrep.2020.05.004

5. Iacopini, I., Petri, G., Barrat, A., & Latora, V. (2019). "Simplicial models of social contagion." Nature Communications 10:2485. DOI: 10.1038/s41467-019-10431-6

6. Petri, G., Expert, P., Turkheimer, F., et al. (2014). "Homological scaffolds of brain functional networks." Journal of the Royal Society Interface 11:20140873. DOI: 10.1098/rsif.2014.0873

7. Gambuzza, L. V., Di Patti, F., Gallo, L., et al. (2021). "Stability of synchronization in simplicial complexes." Nature Communications 12:1255. DOI: 10.1038/s41467-021-21486-9

8. Reimann, M. W., Nolte, M., Scolamiero, M., et al. (2017). "Cliques of neurons bound into cavities provide a missing link between structure and function." Frontiers in Computational Neuroscience 11:48. DOI: 10.3389/fncom.2017.00048

9. Millán, A. P., Torres, J. J., & Bianconi, G. (2020). "Explosive higher-order Kuramoto dynamics on simplicial complexes." Physical Review Letters 124:218301. DOI: 10.1103/PhysRevLett.124.218301

10. 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

---

Numerical Considerations

Euler vs Higher-Order Methods

The HypergraphEngine uses Euler integration. For stiff systems (large \(\sigma\) or mixed orders), the time step must satisfy:

\[\Delta t < \frac{1}{\max_i |\dot{\theta}_i|} \approx \frac{1}{\max(\omega) + \sigma_{max} \cdot N}\]

For \(N = 50\), \(\sigma = 1\): \(\Delta t < 0.02\) is typically safe. For higher accuracy, reduce \(\Delta t\) or use SplittingEngine for the pairwise component combined with explicit hyperedge steps.

Floating-Point Precision

The phase sum \(\Phi_e = \sum_{j \in e} \theta_j\) for large \(k\) can accumulate rounding errors. For \(k > 10\), the argument to \(\sin\) may lose 1-2 digits of precision. This is rarely an issue in practice since most physical hyperedges have \(k \leq 5\).

Determinism

Both Python and Rust paths are deterministic given the same inputs. The Rust path uses the same mathematical operations in the same order, ensuring bitwise-identical results across runs.

Direct Mojo hypergraph adapters enforce exact raw stdout cardinality and finite phase-domain values: HGRUN must emit exactly one finite scalar in [0, 2*pi) per oscillator. Blank, truncated, overlong, non-numeric, non-finite, or out-of-domain output is rejected before public arrays are returned.

---

Edge Cases and Limitations

Empty Hyperedge List

With no hyperedges added, run() uses the Python loop (Rust path requires at least one hyperedge). If pairwise knm is provided, the engine acts as a standard Kuramoto integrator (but without RK4 — use UPDEEngine for higher-order integration).

Single-Node Hyperedges

A hyperedge with \(k = 1\) is mathematically trivial: the argument \(\Phi_e - k\theta_m = \theta_m - \theta_m = 0\), so \(\sin(0) = 0\). No contribution. The engine does not filter these, but they have zero effect.

Very Large \(|E|\)

For \(N = 100\) with all-to-all 3-body: \(C(100, 3) = 161{,}700\) edges. Each step requires one \(\sin\) evaluation per node per edge, giving \(\sim 16\) million trig calls per step. At \(N > 50\), consider whether SimplicialEngine (all-to-all implicit at \(O(N^2)\)) is more appropriate.

Duplicate Edges

Adding the same edge twice doubles its effective coupling strength. The engine does not deduplicate — this is by design, as different edges between the same nodes may represent distinct physical interactions.

---

Appendix: Relationship to Other Engines

If you need...	Use...
Only pairwise	UPDEEngine
Only 3-body (all-to-all)	SimplicialEngine (faster)
Any \(k\)-body (explicit edges)	HypergraphEngine (this)
Mixed pairwise + 3-body	Either (Simplicial is faster)
Sparse edge sets	HypergraphEngine
Dense all-to-all 3-body	SimplicialEngine (\(O(N^2)\) vs \(O(N^3)\))

---

Appendix B: Choosing the Right Engine

Decision tree for higher-order Kuramoto integration:

Is the coupling purely pairwise?
  ├─ Yes → UPDEEngine (RK4/RK45, fastest for pairwise)
  └─ No → Are ALL higher-order interactions exactly 3-body?
           ├─ Yes → Is it all-to-all (every triplet)?
           │         ├─ Yes → SimplicialEngine (O(N²) factorisation)
           │         └─ No  → HypergraphEngine (explicit edge list)
           └─ No → HypergraphEngine (mixed k-body, this module)

For systems with \(N > 100\) and dense 3-body coupling, always prefer SimplicialEngine — the \(O(N^2)\) vs \(O(N^3)\) difference is significant (e.g. \(N=100\): 10,000 vs 166,650 operations per step).

---

Appendix C: Complexity Summary

Component	Python	Rust	Notes
Pairwise coupling	\(O(N^2 \cdot T)\)	\(O(N^2 \cdot T)\)	Same as UPDEEngine
Per-edge coupling	$O(	E	\cdot k \cdot T)$
All-to-all \(k=3\)	\(O(N^3 \cdot T / 6)\)	\(O(N^3 \cdot T / 6)\)	Use SimplicialEngine for \(O(N^2)\)
Edge encoding (FFI)	$O(	E	\cdot k)$
Memory (edges)	$O(	E	\cdot k)$
Memory (state)	\(O(N)\)	\(O(N)\)	Phase vector