diff --git a/4-Infrastructure/shim/capability_probe.py b/4-Infrastructure/shim/capability_probe.py
new file mode 100644
index 00000000..7a3ffac7
--- /dev/null
+++ b/4-Infrastructure/shim/capability_probe.py
@@ -0,0 +1,405 @@
+#!/usr/bin/env python3
+"""Capability probe: detect available quantum/optimization/compute backends.
+
+Probes each package, runs a minimal functional test, and emits a structured
+receipt mapping capability slots → available backends + version + status.
+
+Usage:
+    python3 capability_probe.py                          # stdout
+    python3 capability_probe.py -o probe_receipt.json   # file
+    python3 capability_probe.py --json                   # compact JSON
+"""
+
+import importlib
+import json
+import math
+import sys
+import time
+import traceback
+from dataclasses import dataclass, field, asdict
+from typing import Any, Optional
+
+
+# =========================================================================
+# Probe registry
+# =========================================================================
+
+@dataclass
+class ProbeResult:
+    package: str
+    version: str = ""
+    import_ok: bool = False
+    functional_ok: bool = False
+    error: str = ""
+    test_detail: str = ""
+    elapsed_s: float = 0.0
+
+
+PROBES: list[dict] = []
+
+
+def register(category: str, slot: str, package: str,
+             import_name: Optional[str] = None,
+             min_version: Optional[str] = None,
+             test: Optional[callable] = None):
+    PROBES.append(dict(
+        category=category,
+        slot=slot,
+        package=package,
+        import_name=import_name or package,
+        min_version=min_version,
+        test=test,
+    ))
+
+
+def _try_import(modname: str) -> tuple[bool, str]:
+    try:
+        mod = importlib.import_module(modname)
+        ver = getattr(mod, "__version__", "")
+        return True, ver
+    except Exception:
+        return False, ""
+
+
+def _run_test(modname: str, test_fn: callable) -> tuple[bool, str]:
+    try:
+        mod = importlib.import_module(modname)
+        result = test_fn(mod)
+        return True, result
+    except Exception as e:
+        tb = traceback.format_exc()
+        return False, f"{type(e).__name__}: {e}\n{tb[:200]}"
+
+
+# — Quantum backends —
+
+
+def _test_perceval(mod):
+    from perceval import BasicState, Circuit, components, Processor, SLOSBackend
+    import math
+    c = Circuit(2)
+    c.add(0, components.BS(theta=math.pi / 4))
+    backend = SLOSBackend()
+    proc = Processor(backend, c)
+    input_state = BasicState("|0,1>")
+    proc.with_input(input_state)
+    output = proc.probs()
+    results = output["results"]
+    return f"BS probabilities: {dict(results)}"
+
+
+def _test_quimb(mod):
+    import quimb.tensor as qtn
+    import numpy as np
+    t1 = qtn.Tensor(np.array([[1.0, 0.0], [0.0, 1.0]]), inds=["a", "b"])
+    t2 = qtn.Tensor(np.array([[1.0, 0.5], [0.0, 1.0]]), inds=["b", "c"])
+    tn = qtn.TensorNetwork([t1, t2])
+    result = tn.contract()
+    return f"contracted 2 tensors, result shape={result.shape}"
+
+
+def _test_cirq(mod):
+    import cirq
+    q = cirq.LineQubit(0)
+    circuit = cirq.Circuit([cirq.H(q), cirq.measure(q)])
+    sim = cirq.Simulator()
+    result = sim.run(circuit, repetitions=10)
+    return f"simulated {len(result.measurements)} outcomes"
+
+
+def _test_qiskit(mod):
+    from qiskit import QuantumCircuit
+    from qiskit_aer import AerSimulator
+    qc = QuantumCircuit(2, 2)
+    qc.h(0)
+    qc.cx(0, 1)
+    qc.measure([0, 1], [0, 1])
+    sim = AerSimulator()
+    result = sim.run(qc, shots=10).result()
+    counts = result.get_counts()
+    return f"simulated {len(counts)} distinct outcomes"
+
+
+# — Classical optimization —
+
+
+def _test_highspy(mod):
+    import highspy
+    import numpy as np
+    # Minimize -2x0 - x1 subject to x0 + x1 >= 1, 0 <= x0, x1 <= 2
+    # Optimum: x0=2, x1=0 → obj=-4
+    lp = highspy.HighsLp()
+    lp.num_col_ = 2
+    lp.num_row_ = 1
+    lp.sense_ = highspy.ObjSense.kMinimize
+    lp.col_cost_ = np.array([-2.0, -1.0])
+    lp.col_lower_ = np.array([0.0, 0.0])
+    lp.col_upper_ = np.array([2.0, 2.0])
+    lp.row_lower_ = np.array([1.0])
+    lp.row_upper_ = np.array([float("inf")])
+    lp.a_matrix_ = highspy.HighsSparseMatrix()
+    lp.a_matrix_.start_ = np.array([0, 1, 2], dtype=np.int32)
+    lp.a_matrix_.index_ = np.array([0, 0], dtype=np.int32)
+    lp.a_matrix_.value_ = np.array([1.0, 1.0])
+    lp.a_matrix_.format_ = highspy.MatrixFormat.kColwise
+    h = highspy.Highs()
+    h.passModel(lp)
+    h.run()
+    info = h.getInfo()
+    sol = h.getSolution()
+    x = list(sol.col_value) if hasattr(sol, "col_value") else []
+    return f"LP solved, obj={info.objective_function_value:.2f}, x={x}"
+
+
+def _test_scipy_optimize(mod):
+    from scipy.optimize import minimize
+    result = minimize(lambda x: x[0]**2 + x[1]**2, [1.0, 1.0], method="Nelder-Mead")
+    return f"minimized to ({result.x[0]:.4f}, {result.x[1]:.4f})"
+
+
+def _test_networkx(mod):
+    import networkx as nx
+    g = nx.complete_graph(5)
+    path = nx.shortest_path(g, source=0, target=4)
+    return f"K_5 diameter={nx.diameter(g)}, path={path}"
+
+
+# — GPU compute —
+
+
+def _test_wgpu(mod):
+    import wgpu
+    adapter = wgpu.gpu.request_adapter_sync(power_preference="high-performance")
+    device = adapter.request_device_sync()
+    info = device.adapter_info
+    dev_name = getattr(info, "name", str(info))
+    dev_type = getattr(info, "adapter_type", "?")
+    backend_type = getattr(info, "backend_type", "?")
+    return f"adapter={dev_name}, type={dev_type}, backend={backend_type}"
+
+
+# — Utility —
+
+
+def _test_opt_einsum(mod):
+    import opt_einsum as oe
+    import numpy as np
+    a = np.random.rand(2, 3)
+    b = np.random.rand(3, 4)
+    c = oe.contract("ij,jk->ik", a, b)
+    # Compute path
+    path_info = oe.contract_path("ij,jk->ik", a, b)
+    return f"contraction shape={c.shape}, path={path_info[1]}"
+
+
+# =========================================================================
+# Register all probes
+# =========================================================================
+
+# Photonic quantum
+register("quantum", "photonic_slos", "perceval-quandela",
+         import_name="perceval", test=_test_perceval)
+# Bosonic tensor network
+register("quantum", "bosonic_tn", "quimb", test=_test_quimb)
+# QAOA circuit sim (Cirq)
+register("quantum", "qaoa_cirq", "cirq", test=_test_cirq)
+# QAOA circuit sim (Qiskit)
+register("quantum", "qaoa_qiskit", "qiskit", test=_test_qiskit)
+# Classical MIP
+register("optimization", "mip_highs", "highspy", test=_test_highspy)
+# Classical continuous
+register("optimization", "continuous_scipy", "scipy.optimize",
+         import_name="scipy", test=_test_scipy_optimize)
+# Graph algorithms
+register("optimization", "graph_networkx", "networkx", test=_test_networkx)
+# GPU compute
+register("compute", "gpu_wgpu", "wgpu", test=_test_wgpu)
+# Tensor contraction
+register("compute", "tensor_opt_einsum", "opt_einsum", test=_test_opt_einsum)
+
+
+# =========================================================================
+# Probe runner
+# =========================================================================
+
+def run_all_probes() -> dict[str, list[dict]]:
+    results: dict[str, list[dict]] = {}
+    for entry in PROBES:
+        cat = entry["category"]
+        pkg = entry["package"]
+        modname = entry["import_name"]
+        test_fn = entry["test"]
+
+        t0 = time.time()
+        import_ok, version = _try_import(modname)
+        functional_ok = False
+        test_detail = ""
+        error = ""
+
+        if import_ok and test_fn:
+            functional_ok, test_detail = _run_test(modname, test_fn)
+            if not functional_ok:
+                error = test_detail
+                test_detail = ""
+
+        if not import_ok:
+            error = f"import failed for {modname}"
+
+        elapsed = time.time() - t0
+
+        result = {
+            "slot": entry["slot"],
+            "package": pkg,
+            "version": version,
+            "import_ok": import_ok,
+            "functional_ok": functional_ok,
+            "test_detail": test_detail if functional_ok else "",
+            "error": error,
+            "elapsed_s": round(elapsed, 3),
+        }
+        results.setdefault(cat, []).append(result)
+
+    return results
+
+
+# =========================================================================
+# Capability assignment
+# =========================================================================
+
+def assign_capabilities(probe_results: dict[str, list[dict]]) -> dict:
+    avail = {}
+    for cat, entries in probe_results.items():
+        for e in entries:
+            if e["functional_ok"]:
+                avail[e["slot"]] = {
+                    "package": e["package"],
+                    "version": e["version"],
+                }
+    return avail
+
+
+def assign_formulation_modes(avail: dict) -> list[str]:
+    """Map available backends to RRC-relevant formulation modes."""
+    modes = []
+    if "bosonic_tn" in avail:
+        modes.append("bosonic_tensor_network")
+    if "qaoa_cirq" in avail or "qaoa_qiskit" in avail:
+        modes.append("qaoa_variational")
+    if "mip_highs" in avail:
+        modes.append("qubo_classical_mip")
+        modes.append("qap_classical_mip")
+    if "photonic_slos" in avail:
+        modes.append("photonic_slos")
+    if "continuous_scipy" in avail:
+        modes.append("continuous_optimization")
+    if "graph_networkx" in avail:
+        modes.append("spectral_graph_pipeline")
+    if "gpu_wgpu" in avail:
+        modes.append("gpu_vulkan_compute")
+    if "tensor_opt_einsum" in avail:
+        modes.append("contraction_optimization")
+    return modes
+
+
+def assign_directed_formulation_support(avail: dict) -> dict:
+    """What's available for directed (asymmetric) routing formulations."""
+    support = {}
+    if "mip_highs" in avail:
+        support["qap"] = {
+            "available": True,
+            "comment": "HiGHS handles asymmetric QAP natively (directed costs preserved in MIP)",
+            "recommended": True,
+            "solver": "highspy",
+        }
+    else:
+        support["qap"] = {
+            "available": False,
+            "comment": "no MIP solver — use SA/Levy on symmetrized QUBO instead",
+        }
+    if "qaoa_cirq" in avail or "qaoa_qiskit" in avail:
+        support["qaoa"] = {
+            "available": True,
+            "comment": "QAOA on QUBO — asymmetry lost in x_i x_j = x_j x_i symmetrization",
+            "recommended": False,
+            "solver": "cirq" if "qaoa_cirq" in avail else "qiskit",
+        }
+    return support
+
+
+# =========================================================================
+# Receipt assembly
+# =========================================================================
+
+def build_receipt(probe_results: dict[str, list[dict]],
+                  avail: dict, modes: list[str],
+                  directed: dict) -> dict:
+    timestamp = time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime())
+    total_probes = sum(len(v) for v in probe_results.values())
+    functional_count = sum(
+        1 for entries in probe_results.values()
+        for e in entries if e["functional_ok"]
+    )
+
+    return {
+        "schema": "capability_probe_receipt_v1",
+        "generated_at_utc": timestamp,
+        "probe_summary": {
+            "total_probes": total_probes,
+            "functional": functional_count,
+            "failed": total_probes - functional_count,
+        },
+        "results": probe_results,
+        "available": avail,
+        "formulation_modes": modes,
+        "directed_formulation_support": directed,
+        "claim_boundary": "capability-probe-only;no-decision-logic",
+    }
+
+
+# =========================================================================
+# CLI
+# =========================================================================
+
+def main():
+    import argparse
+    parser = argparse.ArgumentParser(description="Capability probe")
+    parser.add_argument("-o", "--output", type=str, help="JSON output path")
+    parser.add_argument("--json", action="store_true", help="Compact JSON (one line)")
+    parser.add_argument("--human", action="store_true", help="Human-readable summary")
+    args = parser.parse_args()
+
+    probe_results = run_all_probes()
+    avail = assign_capabilities(probe_results)
+    modes = assign_formulation_modes(avail)
+    directed = assign_directed_formulation_support(avail)
+    receipt = build_receipt(probe_results, avail, modes, directed)
+
+    if args.output:
+        with open(args.output, "w") as f:
+            json.dump(receipt, f, indent=2)
+        print(f"Wrote receipt to {args.output}")
+
+    if args.json:
+        print(json.dumps(receipt))
+
+    if args.human or not args.json:
+        print(f"\n=== Capability Probe ===")
+        print(f"  {receipt['probe_summary']['functional']}/{receipt['probe_summary']['total_probes']} probes functional")
+        print(f"\n  Available backends ({len(avail)}):")
+        for slot, info in sorted(avail.items()):
+            print(f"    {slot:30s}  {info['package']:20s}  v{info['version']}")
+        print(f"\n  Formulation modes ({len(modes)}):")
+        for m in modes:
+            print(f"    {m}")
+        print(f"\n  Directed (asymmetric) routing support:")
+        for label, info in directed.items():
+            status = "✅" if info["available"] else "❌"
+            rec = " (recommended)" if info.get("recommended") else " (not recommended)" if info.get("recommended") is False else ""
+            print(f"    {status} {label}: {info['comment']}{rec}")
+
+    if not (args.output or args.json or args.human):
+        print(json.dumps(receipt, indent=2))
+
+
+if __name__ == "__main__":
+    main()
diff --git a/shared-data/artifacts/capability_probe_receipt_2026-06-21.json b/shared-data/artifacts/capability_probe_receipt_2026-06-21.json
new file mode 100644
index 00000000..53f6cf7b
--- /dev/null
+++ b/shared-data/artifacts/capability_probe_receipt_2026-06-21.json
@@ -0,0 +1,171 @@
+{
+  "schema": "capability_probe_receipt_v1",
+  "generated_at_utc": "2026-06-21T04:51:09Z",
+  "probe_summary": {
+    "total_probes": 9,
+    "functional": 9,
+    "failed": 0
+  },
+  "results": {
+    "quantum": [
+      {
+        "slot": "photonic_slos",
+        "package": "perceval-quandela",
+        "version": "1.2.3",
+        "import_ok": true,
+        "functional_ok": true,
+        "test_detail": "BS probabilities: {|0,1>: 0.8535533905932737, |1,0>: 0.14644660940672624}",
+        "error": "",
+        "elapsed_s": 0.536
+      },
+      {
+        "slot": "bosonic_tn",
+        "package": "quimb",
+        "version": "1.14.0",
+        "import_ok": true,
+        "functional_ok": true,
+        "test_detail": "contracted 2 tensors, result shape=(2, 2)",
+        "error": "",
+        "elapsed_s": 0.194
+      },
+      {
+        "slot": "qaoa_cirq",
+        "package": "cirq",
+        "version": "1.6.1",
+        "import_ok": true,
+        "functional_ok": true,
+        "test_detail": "simulated 1 outcomes",
+        "error": "",
+        "elapsed_s": 0.231
+      },
+      {
+        "slot": "qaoa_qiskit",
+        "package": "qiskit",
+        "version": "2.4.2",
+        "import_ok": true,
+        "functional_ok": true,
+        "test_detail": "simulated 2 distinct outcomes",
+        "error": "",
+        "elapsed_s": 0.074
+      }
+    ],
+    "optimization": [
+      {
+        "slot": "mip_highs",
+        "package": "highspy",
+        "version": "",
+        "import_ok": true,
+        "functional_ok": true,
+        "test_detail": "LP solved, obj=-6.00, x=[2.0, 2.0]",
+        "error": "",
+        "elapsed_s": 0.006
+      },
+      {
+        "slot": "continuous_scipy",
+        "package": "scipy.optimize",
+        "version": "1.18.0",
+        "import_ok": true,
+        "functional_ok": true,
+        "test_detail": "minimized to (-0.0000, 0.0000)",
+        "error": "",
+        "elapsed_s": 0.001
+      },
+      {
+        "slot": "graph_networkx",
+        "package": "networkx",
+        "version": "3.6.1",
+        "import_ok": true,
+        "functional_ok": true,
+        "test_detail": "K_5 diameter=1, path=[0, 4]",
+        "error": "",
+        "elapsed_s": 0.001
+      }
+    ],
+    "compute": [
+      {
+        "slot": "gpu_wgpu",
+        "package": "wgpu",
+        "version": "0.31.0",
+        "import_ok": true,
+        "functional_ok": true,
+        "test_detail": "adapter=<GPUAdapterInfo dict with vendor='NVIDIA', architecture='', device='NVIDIA GeForce RTX 4070 SUPER', description='610.43.02', stackgroup_min_size=4, stackgroup_max_size=128, vendor_id=4318, device_id=10115, adapter_type='DiscreteGPU', backend_type='Vulkan'>, type=?, backend=?",
+        "error": "",
+        "elapsed_s": 0.375
+      },
+      {
+        "slot": "tensor_opt_einsum",
+        "package": "opt_einsum",
+        "version": "3.4.0",
+        "import_ok": true,
+        "functional_ok": true,
+        "test_detail": "contraction shape=(2, 4), path=  Complete contraction:  ij,jk->ik\n         Naive scaling:  3\n     Optimized scaling:  3\n      Naive FLOP count:  4.800e+1\n  Optimized FLOP count:  4.800e+1\n   Theoretical speedup:  1.000e+0\n  Largest intermediate:  8.000e+0 elements\n--------------------------------------------------------------------------------\nscaling        BLAS                current                             remaining\n--------------------------------------------------------------------------------\n   3           GEMM              jk,ij->ik                                ik->ik",
+        "error": "",
+        "elapsed_s": 0.0
+      }
+    ]
+  },
+  "available": {
+    "photonic_slos": {
+      "package": "perceval-quandela",
+      "version": "1.2.3"
+    },
+    "bosonic_tn": {
+      "package": "quimb",
+      "version": "1.14.0"
+    },
+    "qaoa_cirq": {
+      "package": "cirq",
+      "version": "1.6.1"
+    },
+    "qaoa_qiskit": {
+      "package": "qiskit",
+      "version": "2.4.2"
+    },
+    "mip_highs": {
+      "package": "highspy",
+      "version": ""
+    },
+    "continuous_scipy": {
+      "package": "scipy.optimize",
+      "version": "1.18.0"
+    },
+    "graph_networkx": {
+      "package": "networkx",
+      "version": "3.6.1"
+    },
+    "gpu_wgpu": {
+      "package": "wgpu",
+      "version": "0.31.0"
+    },
+    "tensor_opt_einsum": {
+      "package": "opt_einsum",
+      "version": "3.4.0"
+    }
+  },
+  "formulation_modes": [
+    "bosonic_tensor_network",
+    "qaoa_variational",
+    "qubo_classical_mip",
+    "qap_classical_mip",
+    "photonic_slos",
+    "continuous_optimization",
+    "spectral_graph_pipeline",
+    "gpu_vulkan_compute",
+    "contraction_optimization"
+  ],
+  "directed_formulation_support": {
+    "qap": {
+      "available": true,
+      "comment": "HiGHS handles asymmetric QAP natively (directed costs preserved in MIP)",
+      "recommended": true,
+      "solver": "highspy"
+    },
+    "qaoa": {
+      "available": true,
+      "comment": "QAOA on QUBO \u2014 asymmetry lost in x_i x_j = x_j x_i symmetrization",
+      "recommended": false,
+      "solver": "cirq"
+    }
+  },
+  "claim_boundary": "capability-probe-only;no-decision-logic"
+}
\ No newline at end of file