diff --git a/.github/workflows/release-pyapi.yml b/.github/workflows/release-pyapi.yml
new file mode 100644
index 000000000..eea3a1217
--- /dev/null
+++ b/.github/workflows/release-pyapi.yml
@@ -0,0 +1,179 @@
+name: Deploy Maturin wheels
+
+on:
+  push:
+    tags:
+      - "*"
+  workflow_dispatch:
+
+permissions:
+  contents: read
+
+jobs:
+  set-version:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - uses: actions/setup-python@v5
+        with:
+          python-version: ${{ vars.PYTHON_VERSION_WORKFLOWS }}
+      - name: Install tomlkit
+        run: pip install tomlkit
+      - name: Bump versions
+        run: python crates/bump-versions.py ${{ github.ref_name }}
+      - name: Upload patched Cargo files
+        uses: actions/upload-artifact@v4
+        with:
+          name: patched-cargo
+          path: |
+            Cargo.toml
+            crates/**/Cargo.toml
+
+  linux:
+    runs-on: ${{ matrix.platform.runner }}
+    needs: [set-version]
+    strategy:
+      matrix:
+        platform:
+          - runner: ubuntu-latest
+            target: x86_64
+          - runner: ubuntu-latest
+            target: x86
+          - runner: ubuntu-latest
+            target: aarch64
+    steps:
+      - uses: actions/checkout@v4
+      - name: Download patched Cargo files
+        uses: actions/download-artifact@v4
+        with:
+          name: patched-cargo
+          path: .
+      - name: Install 32-bit glibc headers (x86 only)
+        if: matrix.platform.target == 'x86'
+        run: |
+          sudo apt-get update
+          sudo apt-get install -y gcc-multilib
+      - uses: actions/setup-python@v5
+        with:
+          python-version: ${{ vars.PYTHON_VERSION_WORKFLOWS }}
+      - name: Install Rust
+        uses: dtolnay/rust-toolchain@stable
+      - name: Build wheels
+        uses: PyO3/maturin-action@v1
+        with:
+          target: ${{ matrix.platform.target }}
+          args: --release --out dist --find-interpreter -m crates/ekore_py/Cargo.toml
+          sccache: "true"
+          manylinux: auto
+      - name: Upload wheels
+        uses: actions/upload-artifact@v4
+        with:
+          name: wheels-linux-${{ matrix.platform.target }}
+          path: dist
+
+  windows:
+    runs-on: ${{ matrix.platform.runner }}
+    needs: [set-version]
+    strategy:
+      matrix:
+        platform:
+          - runner: windows-latest
+            target: x64
+          - runner: windows-latest
+            target: x86
+    steps:
+      - uses: actions/checkout@v4
+      - name: Download patched Cargo files
+        uses: actions/download-artifact@v4
+        with:
+          name: patched-cargo
+          path: .
+      - uses: actions/setup-python@v5
+        with:
+          python-version: ${{ vars.PYTHON_VERSION_WORKFLOWS }}
+          architecture: ${{ matrix.platform.target }}
+      - name: Install Rust
+        uses: dtolnay/rust-toolchain@stable
+      - name: Build wheels
+        uses: PyO3/maturin-action@v1
+        with:
+          target: ${{ matrix.platform.target }}
+          args: --release --out dist --find-interpreter -m crates/ekore_py/Cargo.toml
+          sccache: "true"
+      - name: Upload wheels
+        uses: actions/upload-artifact@v4
+        with:
+          name: wheels-windows-${{ matrix.platform.target }}
+          path: dist
+
+  macos:
+    runs-on: ${{ matrix.platform.runner }}
+    needs: [set-version]
+    strategy:
+      matrix:
+        platform:
+          - runner: macos-latest
+            target: x86_64
+          - runner: macos-14
+            target: aarch64
+    steps:
+      - uses: actions/checkout@v4
+      - name: Download patched Cargo files
+        uses: actions/download-artifact@v4
+        with:
+          name: patched-cargo
+          path: .
+      - uses: actions/setup-python@v5
+        with:
+          python-version: ${{ vars.PYTHON_VERSION_WORKFLOWS }}
+      - name: Install Rust
+        uses: dtolnay/rust-toolchain@stable
+      - name: Build wheels
+        uses: PyO3/maturin-action@v1
+        with:
+          target: ${{ matrix.platform.target }}
+          args: --release --out dist --find-interpreter -m crates/ekore_py/Cargo.toml
+          sccache: "true"
+      - name: Upload wheels
+        uses: actions/upload-artifact@v4
+        with:
+          name: wheels-macos-${{ matrix.platform.target }}
+          path: dist
+
+  sdist:
+    runs-on: ubuntu-latest
+    needs: [set-version]
+    steps:
+      - uses: actions/checkout@v4
+      - name: Download patched Cargo files
+        uses: actions/download-artifact@v4
+        with:
+          name: patched-cargo
+          path: .
+      - name: Build sdist
+        uses: PyO3/maturin-action@v1
+        with:
+          command: sdist
+          args: --out dist -m crates/ekore_py/Cargo.toml
+      - name: Upload sdist
+        uses: actions/upload-artifact@v4
+        with:
+          name: wheels-sdist
+          path: dist
+
+  release:
+    name: Release
+    runs-on: ubuntu-latest
+    if: "startsWith(github.ref, 'refs/tags/')"
+    needs: [linux, windows, macos, sdist]
+    steps:
+      - uses: actions/download-artifact@v4
+        with:
+          pattern: wheels-*
+          merge-multiple: true
+          path: dist
+      - name: Publish to PyPI
+        uses: pypa/gh-action-pypi-publish@release/v1
+        with:
+          password: ${{ secrets.PYPI_TOKEN }}
+          skip-existing: true
diff --git a/.github/workflows/unittests-pyapi.yml b/.github/workflows/unittests-pyapi.yml
new file mode 100644
index 000000000..80529d5fe
--- /dev/null
+++ b/.github/workflows/unittests-pyapi.yml
@@ -0,0 +1,25 @@
+name: Python-API unit tests
+
+on: push
+
+jobs:
+  test-pyapi:
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v4
+      - name: Setup Python
+        uses: actions/setup-python@v5
+        with:
+          python-version: ${{ vars.PYTHON_VERSION_WORKFLOWS }}
+      - name: Setup Rust toolchain
+        uses: dtolnay/rust-toolchain@stable
+      - name: Create virtual environment
+        run: python -m venv .venv
+      - name: Add venv to PATH
+        run: echo "$PWD/.venv/bin" >> "$GITHUB_PATH"
+      - name: Install task runner and test deps
+        run: pip install maturin numpy poethepoet pytest
+      - name: Run Python-API tests
+        run: |
+          poe build-pyapi
+          poe ptest
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 73b8d2d0c..e9d012337 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,7 +13,12 @@ We may prefix the items to signal the scope:
 
 Items without prefix refer to a global change.
 
-## [Unreleased](https://github.com/NNPDF/eko/compare/v0.15.4...HEAD)
+## [Unreleased](https://github.com/NNPDF/eko/compare/v0.15.5...HEAD)
+
+### Added
+- rust: `ekore_py` crate, which exposes `ekore` crate as a Python bindings ([#554](https://github.com/NNPDF/eko/pull/554))
+
+## [0.15.5](https://github.com/NNPDF/eko/compare/v0.15.4...v0.15.5) - 2026-07-22
 
 ### Added
 - rust: `ekore_capi` crate, which exposes `ekore` crate as a C-ABI for other languages ([#537](https://github.com/NNPDF/eko/pull/537))
diff --git a/Cargo.lock b/Cargo.lock
index 4ac8fc2b8..9d12b7b3d 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -1,6 +1,6 @@
 # This file is automatically @generated by Cargo.
 # It is not intended for manual editing.
-version = 3
+version = 4
 
 [[package]]
 name = "GSL"
@@ -239,6 +239,15 @@ dependencies = [
  "num",
 ]
 
+[[package]]
+name = "ekore_py"
+version = "0.0.1"
+dependencies = [
+ "ekore",
+ "numpy",
+ "pyo3",
+]
+
 [[package]]
 name = "encoding_rs"
 version = "0.8.35"
@@ -369,6 +378,12 @@ dependencies = [
  "hashbrown",
 ]
 
+[[package]]
+name = "heck"
+version = "0.5.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2304e00983f87ffb38b55b444b5e3b60a884b5d30c0fca7d82fe33449bbe55ea"
+
 [[package]]
 name = "ignore"
 version = "0.4.23"
@@ -554,11 +569,27 @@ dependencies = [
  "autocfg",
 ]
 
+[[package]]
+name = "numpy"
+version = "0.29.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6a5b15d63a5ff39e378daed0e1340d3a5964703ea9712eb09a0dc66fade996f4"
+dependencies = [
+ "libc",
+ "ndarray",
+ "num-complex",
+ "num-integer",
+ "num-traits",
+ "pyo3",
+ "pyo3-build-config",
+ "rustc-hash",
+]
+
 [[package]]
 name = "once_cell"
-version = "1.20.3"
+version = "1.21.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "945462a4b81e43c4e3ba96bd7b49d834c6f61198356aa858733bc4acf3cbe62e"
+checksum = "9f7c3e4beb33f85d45ae3e3a1792185706c8e16d043238c593331cc7cd313b50"
 
 [[package]]
 name = "paste"
@@ -617,6 +648,12 @@ version = "0.3.31"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "953ec861398dccce10c670dfeaf3ec4911ca479e9c02154b3a215178c5f566f2"
 
+[[package]]
+name = "portable-atomic"
+version = "1.14.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "3d20d5497ef88037a52ff98267d066e7f11fcc5e99bbfbd58a42336193aacec3"
+
 [[package]]
 name = "predicates"
 version = "3.1.3"
@@ -669,6 +706,64 @@ dependencies = [
  "pest_derive",
 ]
 
+[[package]]
+name = "pyo3"
+version = "0.29.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "cd274650b21d4bfc26a0a47587962c1edb425f69287324355cd040c3ea66071c"
+dependencies = [
+ "libc",
+ "num-complex",
+ "once_cell",
+ "portable-atomic",
+ "pyo3-build-config",
+ "pyo3-ffi",
+ "pyo3-macros",
+]
+
+[[package]]
+name = "pyo3-build-config"
+version = "0.29.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "c5e2a7d2f0d013342f295c048ad19237add5154a55b1c5a254c0ec93d4109078"
+dependencies = [
+ "target-lexicon",
+]
+
+[[package]]
+name = "pyo3-ffi"
+version = "0.29.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ca85c467da1bbc8d866eea5deff9cf29ea5f7785054a17da36e65bda9c05845b"
+dependencies = [
+ "libc",
+ "pyo3-build-config",
+]
+
+[[package]]
+name = "pyo3-macros"
+version = "0.29.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "9ac53762fd065daa3194dd09337a38bd793a188100fd1a9304c4ab312d901771"
+dependencies = [
+ "proc-macro2",
+ "pyo3-macros-backend",
+ "quote",
+ "syn",
+]
+
+[[package]]
+name = "pyo3-macros-backend"
+version = "0.29.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "4ca3a1557399783172dc5bf39cfca835157732532cba56b71d2292161e53b362"
+dependencies = [
+ "heck",
+ "proc-macro2",
+ "quote",
+ "syn",
+]
+
 [[package]]
 name = "quote"
 version = "1.0.38"
@@ -722,6 +817,12 @@ version = "0.8.5"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "2b15c43186be67a4fd63bee50d0303afffcef381492ebe2c5d87f324e1b8815c"
 
+[[package]]
+name = "rustc-hash"
+version = "2.1.3"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "6b1e7f9a428571be2dc5bc0505c13fb6bf936822b894ec87abf8a08a4e51742d"
+
 [[package]]
 name = "rustix"
 version = "0.38.44"
@@ -797,6 +898,12 @@ dependencies = [
  "xattr",
 ]
 
+[[package]]
+name = "target-lexicon"
+version = "0.13.5"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "adb6935a6f5c20170eeceb1a3835a49e12e19d792f6dd344ccc76a985ca5a6ca"
+
 [[package]]
 name = "tempfile"
 version = "3.17.1"
diff --git a/Cargo.toml b/Cargo.toml
index 39ac3d239..3716c60e7 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -31,6 +31,8 @@ lz4_flex = "0.11.6"
 ndarray = "0.15.4"
 ndarray-npy = "0.8.1"
 num = "0.4.1"
+numpy = "0.29.0"
+pyo3 = "0.29.0"
 tar = "0.4.46"
 thiserror = "1.0.63"
 yaml-rust2 = "0.8"
diff --git a/crates/dekoder/README.md b/crates/dekoder/README.md
index 086b29de8..ed9e34dac 100644
--- a/crates/dekoder/README.md
+++ b/crates/dekoder/README.md
@@ -17,6 +17,7 @@ When using our code please cite
 - [eko](https://crates.io/crates/eko) - Core EKO utilities
 - [ekore](https://crates.io/crates/ekore) - Anomalous dimensions and operator matrix elements
 - [ekore_capi](https://crates.io/crates/ekore_capi) - C API for ekore
+- [ekore_py](https://crates.io/crates/ekore_py) - Python API for ekore
 
 ## License
 
diff --git a/crates/eko/README.md b/crates/eko/README.md
index 22b475517..b884cee6a 100644
--- a/crates/eko/README.md
+++ b/crates/eko/README.md
@@ -17,6 +17,7 @@ When using our code please cite
 - [eko](https://crates.io/crates/eko) - Core EKO utilities
 - [ekore](https://crates.io/crates/ekore) - Anomalous dimensions and operator matrix elements
 - [ekore_capi](https://crates.io/crates/ekore_capi) - C API for ekore
+- [ekore_py](https://crates.io/crates/ekore_py) - Python API for ekore
 
 ## License
 
diff --git a/crates/ekore/README.md b/crates/ekore/README.md
index 4047e378f..87fb6fe92 100644
--- a/crates/ekore/README.md
+++ b/crates/ekore/README.md
@@ -18,6 +18,7 @@ When using our code please cite
 - [eko](https://crates.io/crates/eko) - Core EKO utilities
 - [ekore](https://crates.io/crates/ekore) - Anomalous dimensions and operator matrix elements
 - [ekore_capi](https://crates.io/crates/ekore_capi) - C API for ekore
+- [ekore_py](https://crates.io/crates/ekore_py) - Python API for ekore
 
 ## License
 
diff --git a/crates/ekore_capi/README.md b/crates/ekore_capi/README.md
index 00f33944c..1a6491d1b 100644
--- a/crates/ekore_capi/README.md
+++ b/crates/ekore_capi/README.md
@@ -17,6 +17,7 @@ When using our code please cite
 - [eko](https://crates.io/crates/eko) - Core EKO utilities
 - [ekore](https://crates.io/crates/ekore) - Anomalous dimensions and operator matrix elements
 - [ekore_capi](https://crates.io/crates/ekore_capi) - C API for ekore
+- [ekore_py](https://crates.io/crates/ekore_py) - Python API for ekore
 
 ## License
 
diff --git a/crates/ekore_py/Cargo.toml b/crates/ekore_py/Cargo.toml
new file mode 100644
index 000000000..9f41e0d04
--- /dev/null
+++ b/crates/ekore_py/Cargo.toml
@@ -0,0 +1,26 @@
+[package]
+name = "ekore_py"
+description = "Python API for the ekore crate"
+readme = "README.md"
+
+authors.workspace = true
+categories.workspace = true
+edition.workspace = true
+keywords.workspace = true
+license.workspace = true
+repository.workspace = true
+rust-version.workspace = true
+version.workspace = true
+exclude = ["tests"]
+
+[package.metadata.docs.rs]
+rustdoc-args = ["--html-in-header", "doc-header.html"]
+
+[lib]
+name = "ekore_py"
+crate-type = ["cdylib"]
+
+[dependencies]
+ekore.workspace = true
+numpy.workspace = true
+pyo3 = { workspace = true, features = ["extension-module", "num-complex"] }
diff --git a/crates/ekore_py/README.md b/crates/ekore_py/README.md
new file mode 100644
index 000000000..1c5499d48
--- /dev/null
+++ b/crates/ekore_py/README.md
@@ -0,0 +1,24 @@
+# ekore_py
+
+Python API for the [ekore](https://crates.io/crates/ekore) crate which is a part of the [EKO](https://github.com/NNPDF/eko) framework.
+
+This crate exposes the anomalous dimensions and operator matrix elements from `ekore` through PyO3, making them callable from Python. The Python extension module is built and packaged using [maturin](https://github.com/PyO3/maturin).
+
+## Citation policy
+
+When using our code please cite
+
+- our DOI: [![DOI](https://zenodo.org/badge/DOI/10.5281/zenodo.3874237.svg)](https://doi.org/10.5281/zenodo.3874237)
+- our paper: [![arXiv](https://img.shields.io/badge/arXiv-2202.02338-b31b1b?labelColor=222222)](https://arxiv.org/abs/2202.02338)
+
+## Crates in the eko framework
+
+- [dekoder](https://crates.io/crates/dekoder) - Reading and writing EKO output files
+- [eko](https://crates.io/crates/eko) - Core EKO utilities
+- [ekore](https://crates.io/crates/ekore) - Anomalous dimensions and operator matrix elements
+- [ekore_capi](https://crates.io/crates/ekore_capi) - C API for ekore
+- [ekore_py](https://crates.io/crates/ekore_py) - Python API for ekore
+
+## License
+
+[GPL-3.0-or-later](https://github.com/NNPDF/eko/blob/master/LICENSE)
diff --git a/crates/ekore_py/doc-header.html b/crates/ekore_py/doc-header.html
new file mode 120000
index 000000000..22282661d
--- /dev/null
+++ b/crates/ekore_py/doc-header.html
@@ -0,0 +1 @@
+../doc-header.html
\ No newline at end of file
diff --git a/crates/ekore_py/pyproject.toml b/crates/ekore_py/pyproject.toml
new file mode 100644
index 000000000..6b270e683
--- /dev/null
+++ b/crates/ekore_py/pyproject.toml
@@ -0,0 +1,17 @@
+[build-system]
+requires = ["maturin>=1.1,<2.0"]
+build-backend = "maturin"
+
+[project]
+name = "ekore_py"
+dynamic = ["version"]
+requires-python = ">=3.11"
+classifiers = [
+  "Programming Language :: Rust",
+  "Programming Language :: Python :: Implementation :: CPython",
+  "Programming Language :: Python :: Implementation :: PyPy",
+]
+
+[tool.pytest.ini_options]
+testpaths = ["tests"]
+addopts = ["-v"]
diff --git a/crates/ekore_py/src/ad_ps.rs b/crates/ekore_py/src/ad_ps.rs
new file mode 100644
index 000000000..96031aa57
--- /dev/null
+++ b/crates/ekore_py/src/ad_ps.rs
@@ -0,0 +1,73 @@
+//! The polarized, space-like anomalous dimensions.
+
+use ekore::anomalous_dimensions::polarized::spacelike;
+use numpy::{Complex64, PyArray1, PyArray3, PyArrayMethods};
+use pyo3::exceptions::PyValueError;
+use pyo3::prelude::*;
+use pyo3::wrap_pyfunction;
+
+use crate::cache::Cache;
+
+/// Compute the tower of non-singlet |QCD| anomalous dimensions.
+///
+/// # Parameters
+/// * `order_qcd`: The QCD coupling power (supported range: < 3).
+/// * `mode`: The specific non-singlet sector.
+/// * `cache`: Harmonic sums cache.
+/// * `nf`: Number of active flavors.
+///
+/// # Returns
+/// A `numpy.ndarray`, complex array of shape `(order_qcd,)`.
+#[pyfunction]
+#[pyo3(signature = (order_qcd, mode, cache, nf))]
+pub fn gamma_ns_qcd<'py>(
+    py: Python<'py>,
+    order_qcd: usize,
+    mode: u16,
+    cache: &Bound<'py, Cache>,
+    nf: u8,
+) -> PyResult<Bound<'py, PyArray1<Complex64>>> {
+    gamma_ns_qcd_body!(
+        py,
+        order_qcd,
+        mode,
+        cache,
+        nf,
+        order_qcd >= 3,
+        spacelike::gamma_ns_qcd
+    )
+}
+
+/// Compute the tower of singlet |QCD| anomalous dimension matrices.
+///
+/// # Parameters
+/// * `order_qcd`: The QCD coupling power (supported range: < 3).
+/// * `cache`: Harmonic sums cache.
+/// * `nf`: Number of active flavors.
+///
+/// # Returns
+/// A `numpy.ndarray`, complex array of shape `(order_qcd, 2, 2)`.
+#[pyfunction]
+#[pyo3(signature = (order_qcd, cache, nf))]
+pub fn gamma_singlet_qcd<'py>(
+    py: Python<'py>,
+    order_qcd: usize,
+    cache: &Bound<'py, Cache>,
+    nf: u8,
+) -> PyResult<Bound<'py, PyArray3<Complex64>>> {
+    gamma_singlet_qcd_body!(
+        py,
+        order_qcd,
+        cache,
+        nf,
+        order_qcd >= 3,
+        spacelike::gamma_singlet_qcd,
+        2
+    )
+}
+
+pub(super) fn register(m: &Bound<'_, PyModule>) -> PyResult<()> {
+    m.add_function(wrap_pyfunction!(gamma_ns_qcd, m)?)?;
+    m.add_function(wrap_pyfunction!(gamma_singlet_qcd, m)?)?;
+    Ok(())
+}
diff --git a/crates/ekore_py/src/ad_us.rs b/crates/ekore_py/src/ad_us.rs
new file mode 100644
index 000000000..077d844d7
--- /dev/null
+++ b/crates/ekore_py/src/ad_us.rs
@@ -0,0 +1,187 @@
+//! The unpolarized, space-like anomalous dimensions.
+
+use ekore::anomalous_dimensions::unpolarized::spacelike;
+use ekore::constants::{MAX_ORDER_QCD, MAX_ORDER_QED};
+use numpy::{Complex64, PyArray1, PyArray2, PyArray3, PyArray4, PyArrayMethods};
+use pyo3::exceptions::PyValueError;
+use pyo3::prelude::*;
+use pyo3::wrap_pyfunction;
+
+use crate::cache::Cache;
+
+/// Compute the tower of the non-singlet |QCD| anomalous dimensions.
+///
+/// # Parameters
+/// * `order_qcd`: The QCD coupling power (must be `<= MAX_ORDER_QCD`).
+/// * `mode`: The specific non-singlet sector.
+/// * `cache`: Harmonic sums cache.
+/// * `nf`: Number of active flavors.
+/// * `n3lo_variation`: The three N3LO variation flags.
+///
+/// # Returns
+/// A `numpy.ndarray`, complex array of shape `(order_qcd,)`.
+#[pyfunction]
+#[pyo3(signature = (order_qcd, mode, cache, nf, n3lo_variation))]
+pub fn gamma_ns_qcd<'py>(
+    py: Python<'py>,
+    order_qcd: usize,
+    mode: u16,
+    cache: &Bound<'py, Cache>,
+    nf: u8,
+    n3lo_variation: [u8; 3],
+) -> PyResult<Bound<'py, PyArray1<Complex64>>> {
+    gamma_ns_qcd_body!(
+        py,
+        order_qcd,
+        mode,
+        cache,
+        nf,
+        n3lo_variation,
+        order_qcd > MAX_ORDER_QCD,
+        spacelike::gamma_ns_qcd
+    )
+}
+
+/// Compute the tower of the singlet |QCD| anomalous dimension matrices.
+///
+/// # Parameters
+/// * `order_qcd`: The QCD coupling power (must be `<= MAX_ORDER_QCD`).
+/// * `cache`: Harmonic sums cache.
+/// * `nf`: Number of active flavors.
+/// * `n3lo_variation`: The four N3LO variation flags.
+///
+/// # Returns
+/// A `numpy.ndarray`, complex array of shape `(order_qcd, 2, 2)`.
+#[pyfunction]
+#[pyo3(signature = (order_qcd, cache, nf, n3lo_variation))]
+pub fn gamma_singlet_qcd<'py>(
+    py: Python<'py>,
+    order_qcd: usize,
+    cache: &Bound<'py, Cache>,
+    nf: u8,
+    n3lo_variation: [u8; 4],
+) -> PyResult<Bound<'py, PyArray3<Complex64>>> {
+    gamma_singlet_qcd_body!(
+        py,
+        order_qcd,
+        cache,
+        nf,
+        n3lo_variation,
+        order_qcd > MAX_ORDER_QCD,
+        spacelike::gamma_singlet_qcd,
+        2
+    )
+}
+
+/// Compute the tower of the |QCD| x |QED| non-singlet anomalous dimensions.
+///
+/// # Parameters
+/// * `order_qcd`: The QCD coupling power (must be `<= MAX_ORDER_QCD`).
+/// * `order_qed`: The QED coupling power (must be `<= MAX_ORDER_QED`).
+/// * `mode`: The specific non-singlet sector.
+/// * `cache`: Harmonic sums cache.
+/// * `nf`: Number of active flavors.
+/// * `n3lo_variation`: The three N3LO variation flags.
+///
+/// # Returns
+/// A `numpy.ndarray`, complex array of shape `(order_qcd + 1, order_qed + 1)`.
+#[pyfunction]
+#[pyo3(signature = (order_qcd, order_qed, mode, cache, nf, n3lo_variation))]
+pub fn gamma_ns_qed<'py>(
+    py: Python<'py>,
+    order_qcd: usize,
+    order_qed: usize,
+    mode: u16,
+    cache: &Bound<'py, Cache>,
+    nf: u8,
+    n3lo_variation: [u8; 3],
+) -> PyResult<Bound<'py, PyArray2<Complex64>>> {
+    gamma_ns_qed_body!(
+        py,
+        order_qcd,
+        order_qed,
+        mode,
+        cache,
+        nf,
+        n3lo_variation,
+        order_qcd > MAX_ORDER_QCD || order_qed > MAX_ORDER_QED,
+        spacelike::gamma_ns_qed
+    )
+}
+
+/// Compute the tower of the |QCD| x |QED| singlet anomalous dimension matrices.
+///
+/// # Parameters
+/// * `order_qcd`: The QCD coupling power (must be `<= MAX_ORDER_QCD`).
+/// * `order_qed`: The QED coupling power (must be `<= MAX_ORDER_QED`).
+/// * `cache`: Harmonic sums cache.
+/// * `nf`: Number of active flavors.
+/// * `n3lo_variation`: The seven N3LO variation flags.
+///
+/// # Returns
+/// A `numpy.ndarray`, complex array of shape `(order_qcd + 1, order_qed + 1, 4, 4)`.
+#[pyfunction]
+#[pyo3(signature = (order_qcd, order_qed, cache, nf, n3lo_variation))]
+pub fn gamma_singlet_qed<'py>(
+    py: Python<'py>,
+    order_qcd: usize,
+    order_qed: usize,
+    cache: &Bound<'py, Cache>,
+    nf: u8,
+    n3lo_variation: [u8; 7],
+) -> PyResult<Bound<'py, PyArray4<Complex64>>> {
+    gamma_qed_matrix_body!(
+        py,
+        order_qcd,
+        order_qed,
+        cache,
+        nf,
+        n3lo_variation,
+        order_qcd > MAX_ORDER_QCD || order_qed > MAX_ORDER_QED,
+        spacelike::gamma_singlet_qed,
+        4
+    )
+}
+
+/// Compute the tower of the |QCD| x |QED| valence anomalous dimension matrices.
+///
+/// # Parameters
+/// * `order_qcd`: The QCD coupling power (must be `<= MAX_ORDER_QCD`).
+/// * `order_qed`: The QED coupling power (must be `<= MAX_ORDER_QED`).
+/// * `cache`: Harmonic sums cache.
+/// * `nf`: Number of active flavors.
+/// * `n3lo_variation`: The three N3LO variation flags.
+///
+/// # Returns
+/// A `numpy.ndarray`, complex array of shape `(order_qcd + 1, order_qed + 1, 2, 2)`.
+#[pyfunction]
+#[pyo3(signature = (order_qcd, order_qed, cache, nf, n3lo_variation))]
+pub fn gamma_valence_qed<'py>(
+    py: Python<'py>,
+    order_qcd: usize,
+    order_qed: usize,
+    cache: &Bound<'py, Cache>,
+    nf: u8,
+    n3lo_variation: [u8; 3],
+) -> PyResult<Bound<'py, PyArray4<Complex64>>> {
+    gamma_qed_matrix_body!(
+        py,
+        order_qcd,
+        order_qed,
+        cache,
+        nf,
+        n3lo_variation,
+        order_qcd > MAX_ORDER_QCD || order_qed > MAX_ORDER_QED,
+        spacelike::gamma_valence_qed,
+        2
+    )
+}
+
+pub(super) fn register(m: &Bound<'_, PyModule>) -> PyResult<()> {
+    m.add_function(wrap_pyfunction!(gamma_ns_qcd, m)?)?;
+    m.add_function(wrap_pyfunction!(gamma_singlet_qcd, m)?)?;
+    m.add_function(wrap_pyfunction!(gamma_ns_qed, m)?)?;
+    m.add_function(wrap_pyfunction!(gamma_singlet_qed, m)?)?;
+    m.add_function(wrap_pyfunction!(gamma_valence_qed, m)?)?;
+    Ok(())
+}
diff --git a/crates/ekore_py/src/cache.rs b/crates/ekore_py/src/cache.rs
new file mode 100644
index 000000000..3e453d896
--- /dev/null
+++ b/crates/ekore_py/src/cache.rs
@@ -0,0 +1,40 @@
+//! Python bindings for the Mellin-space harmonics cache.
+
+use ekore::harmonics::cache::Cache as EkoreCache;
+use numpy::Complex64;
+use pyo3::prelude::*;
+
+/// Mellin-space harmonics cache.
+///
+/// Memoizes harmonic sums so repeated evaluations at the same Mellin N don't recompute them.
+///
+/// # Parameters
+/// * `n` : Mellin variable N.
+///
+/// # Returns
+/// * Returns a `Cache` at the given parameter.
+#[pyclass(name = "Cache", module = "ekore_py")]
+pub struct Cache {
+    pub(crate) inner: EkoreCache,
+    n: Complex64,
+}
+
+#[pymethods]
+impl Cache {
+    #[new]
+    fn new(n: Complex64) -> Self {
+        Self {
+            inner: EkoreCache::new(n),
+            n,
+        }
+    }
+
+    #[getter]
+    fn n(&self) -> Complex64 {
+        self.n
+    }
+
+    fn __repr__(&self) -> String {
+        format!("Cache(n=({}{:+}j))", self.n.re, self.n.im)
+    }
+}
diff --git a/crates/ekore_py/src/constants.rs b/crates/ekore_py/src/constants.rs
new file mode 100644
index 000000000..6cd93aa08
--- /dev/null
+++ b/crates/ekore_py/src/constants.rs
@@ -0,0 +1,36 @@
+//! Global constants.
+
+use pyo3::prelude::*;
+
+/// Maximum QCD coupling power implemented.
+pub const MAX_ORDER_QCD: usize = ekore::constants::MAX_ORDER_QCD;
+/// Maximum QED coupling power implemented.
+pub const MAX_ORDER_QED: usize = ekore::constants::MAX_ORDER_QED;
+
+/// singlet-like non-singlet |PID|.
+pub const PID_NSP: u16 = ekore::constants::PID_NSP;
+/// valence-like non-singlet |PID|.
+pub const PID_NSM: u16 = ekore::constants::PID_NSM;
+/// non-singlet all-valence |PID|.
+pub const PID_NSV: u16 = ekore::constants::PID_NSV;
+/// singlet-like non-singlet up-sector |PID|.
+pub const PID_NSP_U: u16 = ekore::constants::PID_NSP_U;
+/// singlet-like non-singlet down-sector |PID|.
+pub const PID_NSP_D: u16 = ekore::constants::PID_NSP_D;
+/// valence-like non-singlet up-sector |PID|.
+pub const PID_NSM_U: u16 = ekore::constants::PID_NSM_U;
+/// valence-like non-singlet down-sector |PID|.
+pub const PID_NSM_D: u16 = ekore::constants::PID_NSM_D;
+
+pub(super) fn register(m: &Bound<'_, PyModule>) -> PyResult<()> {
+    m.add("MAX_ORDER_QCD", MAX_ORDER_QCD)?;
+    m.add("MAX_ORDER_QED", MAX_ORDER_QED)?;
+    m.add("PID_NSP", PID_NSP)?;
+    m.add("PID_NSM", PID_NSM)?;
+    m.add("PID_NSV", PID_NSV)?;
+    m.add("PID_NSP_U", PID_NSP_U)?;
+    m.add("PID_NSP_D", PID_NSP_D)?;
+    m.add("PID_NSM_U", PID_NSM_U)?;
+    m.add("PID_NSM_D", PID_NSM_D)?;
+    Ok(())
+}
diff --git a/crates/ekore_py/src/lib.rs b/crates/ekore_py/src/lib.rs
new file mode 100644
index 000000000..d56ce5a75
--- /dev/null
+++ b/crates/ekore_py/src/lib.rs
@@ -0,0 +1,108 @@
+//! Python bindings for [`ekore`], the crate providing the anomalous dimensions and operator
+//! matrix elements of the [EKO](https://github.com/NNPDF/eko) framework.
+//!
+//! This crate re-exposes those quantities as a native Python extension module, built with
+//! [PyO3](https://pyo3.rs) and packaged with [maturin](https://www.maturin.rs). See the main
+//! [EKO documentation](https://eko.readthedocs.io/en/latest/) for the physics behind the
+//! computed quantities, and the [ekore docs](https://docs.rs/ekore/latest/ekore) for the
+//! underlying Rust API.
+//!
+//! # Building & consuming
+//!
+//! ```sh
+//! pip install maturin
+//! maturin develop --release --manifest-path crates/ekore_py/Cargo.toml
+//! ```
+//!
+//! ```python
+//! import ekore_py
+//!
+//! cache = ekore_py.Cache(2.0 + 0.0j)
+//! print(ekore_py.ad_ps.gamma_ns_qcd(2, ekore_py.constants.PID_NSP, cache, nf=4))
+//! ```
+//!
+//! # Naming convention
+//!
+//! Every function lives in a submodule named `<family>_<sector>`:
+//!
+//! * `ad_us` - **a**nomalous **d**imensions, **u**npolarized, **s**pace-like
+//! * `ad_ps` - **a**nomalous **d**imensions, **p**olarized, **s**pace-like
+//! * `ome_us` - **o**perator **m**atrix **e**lements, **u**npolarized, **s**pace-like
+//!
+//! e.g. `ekore_py.ad_us.gamma_ns_qcd` is the non-singlet |QCD| anomalous dimension from the
+//! unpolarized, space-like sector.
+//!
+//! # Result shape convention
+//!
+//! Quantities depending on a perturbative order are returned as a NumPy array of complex128:
+//! one entry per order along the leading axis(es), matrix-valued quantities carry their
+//! `(dim, dim)` on the trailing axes.
+//!
+//! # Available perturbative orders
+//!
+//! For the list of available perturbative orders and their associated references check the
+//! [ekore docs](https://docs.rs/ekore/latest/ekore).
+
+#[macro_use]
+mod macros;
+
+pub mod ad_ps;
+pub mod ad_us;
+pub mod cache;
+pub mod constants;
+pub mod ome_us;
+
+use pyo3::prelude::*;
+
+/// A Python module implemented in Rust.
+#[pymodule]
+mod ekore_py {
+    use pyo3::prelude::*;
+    use pyo3::types::PyModule;
+
+    #[pymodule_init]
+    fn init(m: &Bound<'_, PyModule>) -> PyResult<()> {
+        let py = m.py();
+
+        // top-level classes
+        m.add_class::<crate::cache::Cache>()?;
+
+        // constants submodule
+        let constants_mod = PyModule::new(py, "constants")?;
+        crate::constants::register(&constants_mod)?;
+        m.add_submodule(&constants_mod)?;
+        // for `from ekore_py.constants import X` to work:
+        py.import("sys")?
+            .getattr("modules")?
+            .set_item("ekore_py.constants", &constants_mod)?;
+
+        // ad_ps submodule
+        let ad_ps_mod = PyModule::new(py, "ad_ps")?;
+        crate::ad_ps::register(&ad_ps_mod)?;
+        m.add_submodule(&ad_ps_mod)?;
+        // for `from ekore_py.ad_ps import X` to work:
+        py.import("sys")?
+            .getattr("modules")?
+            .set_item("ekore_py.ad_ps", &ad_ps_mod)?;
+
+        // ad_us submodule
+        let ad_us_mod = PyModule::new(py, "ad_us")?;
+        crate::ad_us::register(&ad_us_mod)?;
+        m.add_submodule(&ad_us_mod)?;
+        // for `from ekore_py.ad_us import X` to work:
+        py.import("sys")?
+            .getattr("modules")?
+            .set_item("ekore_py.ad_us", &ad_us_mod)?;
+
+        // ome_us submodule
+        let ome_us_mod = PyModule::new(py, "ome_us")?;
+        crate::ome_us::register(&ome_us_mod)?;
+        m.add_submodule(&ome_us_mod)?;
+        // for `from ekore_py.ome_us import X` to work:
+        py.import("sys")?
+            .getattr("modules")?
+            .set_item("ekore_py.ome_us", &ome_us_mod)?;
+
+        Ok(())
+    }
+}
diff --git a/crates/ekore_py/src/macros.rs b/crates/ekore_py/src/macros.rs
new file mode 100644
index 000000000..b966bf8d0
--- /dev/null
+++ b/crates/ekore_py/src/macros.rs
@@ -0,0 +1,216 @@
+//! Shared bodies for the `#[pyfunction]`s defined throughout this crate.
+
+/// Body for a non-singlet tower function returning shape `(order_qcd,)`.
+macro_rules! gamma_ns_qcd_body {
+    // With an `n3lo_variation` argument.
+    (
+        $py:expr, $order_qcd:expr, $mode:expr, $cache:expr, $nf:expr, $n3lo_variation:expr,
+        $bound:expr, $path:path
+    ) => {{
+        gamma_ns_qcd_body!(@check $order_qcd, $mode, $bound);
+        let mut cache = $cache.borrow_mut();
+        let gamma = $path($order_qcd, $mode, &mut cache.inner, $nf, $n3lo_variation);
+        gamma_ns_qcd_body!(@collect $py, $order_qcd, gamma)
+    }};
+    // Without an `n3lo_variation` argument.
+    (
+        $py:expr, $order_qcd:expr, $mode:expr, $cache:expr, $nf:expr,
+        $bound:expr, $path:path
+    ) => {{
+        gamma_ns_qcd_body!(@check $order_qcd, $mode, $bound);
+        let mut cache = $cache.borrow_mut();
+        let gamma = $path($order_qcd, $mode, &mut cache.inner, $nf);
+        gamma_ns_qcd_body!(@collect $py, $order_qcd, gamma)
+    }};
+    (@check $order_qcd:expr, $mode:expr, $bound:expr) => {
+        if $bound {
+            return Err(PyValueError::new_err(format!(
+                "order_qcd out of the supported range, got {}",
+                $order_qcd
+            )));
+        }
+        if !matches!($mode, $crate::constants::PID_NSP | $crate::constants::PID_NSM | $crate::constants::PID_NSV) {
+            return Err(PyValueError::new_err(format!(
+                "invalid non-singlet mode: {}",
+                $mode
+            )));
+        }
+    };
+    (@collect $py:expr, $order_qcd:expr, $gamma:expr) => {{
+        let data: Vec<Complex64> = $gamma
+            .into_iter()
+            .take($order_qcd)
+            .map(|c| Complex64::new(c.re, c.im))
+            .collect();
+        Ok(PyArray1::from_vec($py, data))
+    }};
+}
+
+/// Body for a singlet/valence |QCD| matrix tower function returning shape `(order_qcd, dim, dim)`.
+macro_rules! gamma_singlet_qcd_body {
+    // With an `n3lo_variation` argument.
+    (
+        $py:expr, $order_qcd:expr, $cache:expr, $nf:expr, $n3lo_variation:expr,
+        $bound:expr, $path:path, $dim:expr
+    ) => {{
+        gamma_singlet_qcd_body!(@check $order_qcd, $bound);
+        let mut cache = $cache.borrow_mut();
+        let gamma = $path($order_qcd, &mut cache.inner, $nf, $n3lo_variation);
+        gamma_singlet_qcd_body!(@collect $py, $order_qcd, gamma, $dim)
+    }};
+    // Without an `n3lo_variation` argument.
+    (
+        $py:expr, $order_qcd:expr, $cache:expr, $nf:expr,
+        $bound:expr, $path:path, $dim:expr
+    ) => {{
+        gamma_singlet_qcd_body!(@check $order_qcd, $bound);
+        let mut cache = $cache.borrow_mut();
+        let gamma = $path($order_qcd, &mut cache.inner, $nf);
+        gamma_singlet_qcd_body!(@collect $py, $order_qcd, gamma, $dim)
+    }};
+    (@check $order_qcd:expr, $bound:expr) => {
+        if $bound {
+            return Err(PyValueError::new_err(format!(
+                "order_qcd out of the supported range, got {}",
+                $order_qcd
+            )));
+        }
+    };
+    (@collect $py:expr, $order_qcd:expr, $gamma:expr, $dim:expr) => {{
+        let mut data: Vec<Complex64> = Vec::with_capacity($order_qcd * $dim * $dim);
+        for mat in $gamma.into_iter().take($order_qcd) {
+            for row in mat.iter() {
+                for v in row.iter() {
+                    data.push(Complex64::new(v.re, v.im));
+                }
+            }
+        }
+        PyArray1::from_vec($py, data)
+            .reshape([$order_qcd, $dim, $dim])
+            .map_err(|e| PyValueError::new_err(e.to_string()))
+    }};
+}
+
+/// Body for a |QCD| x |QED| non-singlet tower function returning shape
+/// `(order_qcd + 1, order_qed + 1)`.
+macro_rules! gamma_ns_qed_body {
+    (
+        $py:expr, $order_qcd:expr, $order_qed:expr, $mode:expr, $cache:expr, $nf:expr,
+        $n3lo_variation:expr, $bound:expr, $path:path
+    ) => {{
+        if $bound {
+            return Err(PyValueError::new_err(format!(
+                "order_qcd/order_qed out of the supported range, got {}, {}",
+                $order_qcd, $order_qed
+            )));
+        }
+        if !matches!(
+            $mode,
+            $crate::constants::PID_NSP_U
+                | $crate::constants::PID_NSP_D
+                | $crate::constants::PID_NSM_U
+                | $crate::constants::PID_NSM_D
+                | $crate::constants::PID_NSP
+                | $crate::constants::PID_NSM
+                | $crate::constants::PID_NSV,
+        ) {
+            return Err(PyValueError::new_err(format!(
+                "invalid non-singlet mode: {}",
+                $mode
+            )));
+        }
+
+        let mut cache = $cache.borrow_mut();
+        let gamma = $path(
+            $order_qcd,
+            $order_qed,
+            $mode,
+            &mut cache.inner,
+            $nf,
+            $n3lo_variation,
+        );
+
+        let mut data: Vec<Complex64> = Vec::with_capacity(($order_qcd + 1) * ($order_qed + 1));
+        for row in gamma.into_iter().take($order_qcd + 1) {
+            for v in row.into_iter().take($order_qed + 1) {
+                data.push(Complex64::new(v.re, v.im));
+            }
+        }
+
+        PyArray1::from_vec($py, data)
+            .reshape([$order_qcd + 1, $order_qed + 1])
+            .map_err(|e| PyValueError::new_err(e.to_string()))
+    }};
+}
+
+/// Body for a |QCD| x |QED| matrix tower function (singlet/valence) returning shape
+/// `(order_qcd + 1, order_qed + 1, dim, dim)`.
+macro_rules! gamma_qed_matrix_body {
+    (
+        $py:expr, $order_qcd:expr, $order_qed:expr, $cache:expr, $nf:expr, $n3lo_variation:expr,
+        $bound:expr, $path:path, $dim:expr
+    ) => {{
+        if $bound {
+            return Err(PyValueError::new_err(format!(
+                "order_qcd/order_qed out of the supported range, got {}, {}",
+                $order_qcd, $order_qed
+            )));
+        }
+
+        let mut cache = $cache.borrow_mut();
+        let gamma = $path(
+            $order_qcd,
+            $order_qed,
+            &mut cache.inner,
+            $nf,
+            $n3lo_variation,
+        );
+
+        let mut data: Vec<Complex64> =
+            Vec::with_capacity(($order_qcd + 1) * ($order_qed + 1) * $dim * $dim);
+        for row in gamma.into_iter().take($order_qcd + 1) {
+            for mat in row.into_iter().take($order_qed + 1) {
+                for r in mat.iter() {
+                    for v in r.iter() {
+                        data.push(Complex64::new(v.re, v.im));
+                    }
+                }
+            }
+        }
+
+        PyArray1::from_vec($py, data)
+            .reshape([$order_qcd + 1, $order_qed + 1, $dim, $dim])
+            .map_err(|e| PyValueError::new_err(e.to_string()))
+    }};
+}
+
+/// Body for an |OME| matrix tower function returning shape `(matching_order_qcd, dim, dim)`.
+macro_rules! ome_matrix_body {
+    (
+        $py:expr, $order:expr, $cache:expr, $nf:expr, $l:expr,
+        $bound:expr, $path:path, $dim:expr
+    ) => {{
+        if $bound {
+            return Err(PyValueError::new_err(format!(
+                "matching_order_qcd out of the supported range, got {}",
+                $order
+            )));
+        }
+
+        let mut cache = $cache.borrow_mut();
+        let ome = $path($order, &mut cache.inner, $nf, $l);
+
+        let mut data: Vec<Complex64> = Vec::with_capacity($order * $dim * $dim);
+        for mat in ome.into_iter().take($order) {
+            for row in mat.iter() {
+                for v in row.iter() {
+                    data.push(Complex64::new(v.re, v.im));
+                }
+            }
+        }
+
+        PyArray1::from_vec($py, data)
+            .reshape([$order, $dim, $dim])
+            .map_err(|e| PyValueError::new_err(e.to_string()))
+    }};
+}
diff --git a/crates/ekore_py/src/ome_us.rs b/crates/ekore_py/src/ome_us.rs
new file mode 100644
index 000000000..58c02b33c
--- /dev/null
+++ b/crates/ekore_py/src/ome_us.rs
@@ -0,0 +1,78 @@
+//! The unpolarized, space-like |OME|.
+#![allow(non_snake_case)]
+
+use ekore::operator_matrix_elements::unpolarized::spacelike;
+use numpy::{Complex64, PyArray1, PyArray3, PyArrayMethods};
+use pyo3::exceptions::PyValueError;
+use pyo3::prelude::*;
+use pyo3::wrap_pyfunction;
+
+use crate::cache::Cache;
+
+/// Compute the tower of the singlet |OME|.
+///
+/// # Parameters
+/// * `matching_order_qcd`: The QCD matching order (supported range: < 3).
+/// * `cache`: Harmonic sums cache.
+/// * `nf`: Number of active flavors.
+/// * `L`: The logarithm parameter.
+///
+/// # Returns
+/// A `numpy.ndarray`, complex array of shape `(matching_order_qcd, 3, 3)`.
+#[pyfunction]
+#[pyo3(signature = (matching_order_qcd, cache, nf, L))]
+pub fn A_singlet<'py>(
+    py: Python<'py>,
+    matching_order_qcd: usize,
+    cache: &Bound<'py, Cache>,
+    nf: u8,
+    L: f64,
+) -> PyResult<Bound<'py, PyArray3<Complex64>>> {
+    ome_matrix_body!(
+        py,
+        matching_order_qcd,
+        cache,
+        nf,
+        L,
+        matching_order_qcd >= 3,
+        spacelike::A_singlet,
+        3
+    )
+}
+
+/// Compute the tower of the non-singlet |OME|.
+///
+/// # Parameters
+/// * `matching_order_qcd`: The QCD matching order (supported range: < 3).
+/// * `cache`: Harmonic sums cache.
+/// * `nf`: Number of active flavors.
+/// * `L`: The logarithm parameter.
+///
+/// # Returns
+/// A `numpy.ndarray`, complex array of shape `(matching_order_qcd, 2, 2)`.
+#[pyfunction]
+#[pyo3(signature = (matching_order_qcd, cache, nf, L))]
+pub fn A_non_singlet<'py>(
+    py: Python<'py>,
+    matching_order_qcd: usize,
+    cache: &Bound<'py, Cache>,
+    nf: u8,
+    L: f64,
+) -> PyResult<Bound<'py, PyArray3<Complex64>>> {
+    ome_matrix_body!(
+        py,
+        matching_order_qcd,
+        cache,
+        nf,
+        L,
+        matching_order_qcd >= 3,
+        spacelike::A_non_singlet,
+        2
+    )
+}
+
+pub(super) fn register(m: &Bound<'_, PyModule>) -> PyResult<()> {
+    m.add_function(wrap_pyfunction!(A_singlet, m)?)?;
+    m.add_function(wrap_pyfunction!(A_non_singlet, m)?)?;
+    Ok(())
+}
diff --git a/crates/ekore_py/tests/__init__.py b/crates/ekore_py/tests/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/crates/ekore_py/tests/test_ad_ps.py b/crates/ekore_py/tests/test_ad_ps.py
new file mode 100644
index 000000000..282cb18e7
--- /dev/null
+++ b/crates/ekore_py/tests/test_ad_ps.py
@@ -0,0 +1,76 @@
+import numpy as np
+import pytest
+from ekore_py import Cache, ad_ps
+from ekore_py.constants import PID_NSM_U, PID_NSP
+from numpy.testing import assert_allclose
+
+CF = 4.0 / 3.0
+CA = 3.0
+TR = 1.0 / 2.0
+
+
+@pytest.mark.parametrize("order_qcd, expected_len", [(0, 0), (1, 1), (2, 2)])
+def test_gamma_ns_qcd_shape(order_qcd, expected_len):
+    c = Cache(0.0 + 0.0j)
+    g = ad_ps.gamma_ns_qcd(order_qcd, PID_NSP, c, nf=3)
+    assert g.shape == (expected_len,)
+    assert g.dtype == np.complex128
+
+
+@pytest.mark.parametrize(
+    "order_qcd, expected_shape", [(0, (0, 2, 2)), (1, (1, 2, 2)), (2, (2, 2, 2))]
+)
+def test_gamma_singlet_qcd_shape(order_qcd, expected_shape):
+    c = Cache(0.0 + 0.0j)
+    g = ad_ps.gamma_singlet_qcd(order_qcd, c, nf=3)
+    assert g.shape == expected_shape
+    assert g.dtype == np.complex128
+
+
+def test_order_qcd_out_of_range_raises():
+    c = Cache(0.0 + 0.0j)
+    with pytest.raises(ValueError):
+        ad_ps.gamma_ns_qcd(3, PID_NSP, c, nf=3)
+    with pytest.raises(ValueError):
+        ad_ps.gamma_singlet_qcd(3, c, nf=3)
+
+
+def test_gamma_ns_qcd():
+    nf = 3
+    c = Cache(1.0 + 0.0j)
+
+    r = ad_ps.gamma_ns_qcd(2, PID_NSP, c, nf)
+
+    assert_allclose(r[0], 0.0, atol=1e-14)
+    assert_allclose(r, np.zeros_like(r), atol=2e-6)
+
+
+def test_gamma_singlet_qcd():
+    nf = 5
+    c = Cache(2.0 + 0.0j)
+
+    g = ad_ps.gamma_singlet_qcd(2, c, nf)
+
+    # LO
+    assert_allclose(g[0, 0, 0] + g[0, 1, 0], 4.0 * CF / 3.0, atol=1e-12)
+    assert_allclose(g[0, 0, 1] + g[0, 1, 1], 3.0 + nf / 3.0, atol=1e-12)
+
+    # NLO
+    expected_nlo = (
+        4.0
+        * nf
+        * (
+            0.574074074 * CF
+            - 2.0 * CA * (-7.0 / 18.0 + 1.0 / 6.0 * (5.0 - np.pi**2 / 3.0))
+        )
+        * TR
+    )
+    assert_allclose(-g[1, 0, 1], expected_nlo, atol=1e-9)
+
+
+def test_guards_unknown_mode_raises():
+    nf = 4
+    c = Cache(1.234 + 0.0j)
+
+    with pytest.raises(ValueError):
+        ad_ps.gamma_ns_qcd(2, PID_NSM_U, c, nf)
diff --git a/crates/ekore_py/tests/test_ad_us.py b/crates/ekore_py/tests/test_ad_us.py
new file mode 100644
index 000000000..b60d92e54
--- /dev/null
+++ b/crates/ekore_py/tests/test_ad_us.py
@@ -0,0 +1,221 @@
+import numpy as np
+import pytest
+from ekore_py import Cache, ad_us
+from ekore_py.constants import (
+    PID_NSM,
+    PID_NSM_D,
+    PID_NSM_U,
+    PID_NSP,
+    PID_NSP_D,
+    PID_NSP_U,
+    PID_NSV,
+)
+from numpy.testing import assert_allclose
+
+MAX_ORDER_QCD = 4
+MAX_ORDER_QED = 2
+
+
+@pytest.mark.parametrize(
+    "order_qcd, expected_len", [(0, 0), (1, 1), (2, 2), (3, 3), (4, 4)]
+)
+def test_gamma_ns_qcd_shape(order_qcd, expected_len):
+    c = Cache(0.0 + 0.0j)
+    var = (0, 0, 0)
+    g = ad_us.gamma_ns_qcd(order_qcd, PID_NSP, c, nf=3, n3lo_variation=var)
+    assert g.shape == (expected_len,)
+    assert g.dtype == np.complex128
+
+
+@pytest.mark.parametrize(
+    "order_qcd, expected_shape",
+    [(0, (0, 2, 2)), (1, (1, 2, 2)), (2, (2, 2, 2)), (3, (3, 2, 2)), (4, (4, 2, 2))],
+)
+def test_gamma_singlet_qcd_shape(order_qcd, expected_shape):
+    c = Cache(0.0 + 0.0j)
+    var = (0, 0, 0, 0)
+    g = ad_us.gamma_singlet_qcd(order_qcd, c, nf=3, n3lo_variation=var)
+    assert g.shape == expected_shape
+    assert g.dtype == np.complex128
+
+
+@pytest.mark.parametrize(
+    "order_qcd, order_qed, expected_shape",
+    [(3, 2, (4, 3)), (1, 1, (2, 2))],
+)
+def test_gamma_ns_qed_shape(order_qcd, order_qed, expected_shape):
+    c = Cache(0.0 + 0.0j)
+    var = (0, 0, 0)
+    g = ad_us.gamma_ns_qed(order_qcd, order_qed, PID_NSP_U, c, nf=3, n3lo_variation=var)
+    assert g.shape == expected_shape
+    assert g.dtype == np.complex128
+
+
+@pytest.mark.parametrize(
+    "order_qcd, order_qed, expected_shape",
+    [(3, 2, (4, 3, 4, 4))],
+)
+def test_gamma_singlet_qed_shape(order_qcd, order_qed, expected_shape):
+    c = Cache(0.0 + 0.0j)
+    var = (0, 0, 0, 0, 0, 0, 0)
+    g = ad_us.gamma_singlet_qed(order_qcd, order_qed, c, nf=3, n3lo_variation=var)
+    assert g.shape == expected_shape
+    assert g.dtype == np.complex128
+
+
+@pytest.mark.parametrize(
+    "order_qcd, order_qed, expected_shape",
+    [(3, 2, (4, 3, 2, 2))],
+)
+def test_gamma_valence_qed_shape(order_qcd, order_qed, expected_shape):
+    c = Cache(0.0 + 0.0j)
+    var = (0, 0, 0)
+    g = ad_us.gamma_valence_qed(order_qcd, order_qed, c, nf=3, n3lo_variation=var)
+    assert g.shape == expected_shape
+    assert g.dtype == np.complex128
+
+
+def test_order_out_of_range_raises():
+    c = Cache(0.0 + 0.0j)
+    with pytest.raises(ValueError):
+        ad_us.gamma_ns_qcd(
+            MAX_ORDER_QCD + 1, PID_NSP, c, nf=3, n3lo_variation=(0, 0, 0)
+        )
+    with pytest.raises(ValueError):
+        ad_us.gamma_singlet_qcd(MAX_ORDER_QCD + 1, c, nf=3, n3lo_variation=(0, 0, 0, 0))
+    with pytest.raises(ValueError):
+        ad_us.gamma_ns_qed(
+            MAX_ORDER_QCD + 1, 1, PID_NSP_U, c, nf=4, n3lo_variation=(0, 0, 0)
+        )
+    with pytest.raises(ValueError):
+        ad_us.gamma_ns_qed(
+            1, MAX_ORDER_QED + 1, PID_NSP_U, c, nf=4, n3lo_variation=(0, 0, 0)
+        )
+    with pytest.raises(ValueError):
+        ad_us.gamma_valence_qed(MAX_ORDER_QCD + 1, 1, c, nf=4, n3lo_variation=(0, 0, 0))
+    with pytest.raises(ValueError):
+        ad_us.gamma_singlet_qed(
+            1, MAX_ORDER_QED + 1, c, nf=4, n3lo_variation=(0, 0, 0, 0, 0, 0, 0)
+        )
+
+
+def test_unknown_ns_mode_raises():
+    c = Cache(1.234 + 0.0j)
+    with pytest.raises(ValueError):
+        ad_us.gamma_ns_qed(2, 0, 10106, c, nf=4, n3lo_variation=(0, 0, 0))
+
+
+def test_gamma_ns_qcd():
+    nf3, nf5 = 3, 5
+    var = (0, 0, 0)
+    c = Cache(1.0 + 0.0j)
+    nsm_refs = [0.06776363, 0.064837, 0.07069]
+    nss_refs = [-0.01100459, -0.00779938, -0.0142098]
+
+    r3 = ad_us.gamma_ns_qcd(3, PID_NSP, c, nf3, var)
+    assert_allclose(r3[0], 0.0, atol=1e-14)
+
+    r2 = ad_us.gamma_ns_qcd(2, PID_NSM, c, nf3, var)
+    assert_allclose(r2, np.zeros_like(r2), atol=2e-6)
+
+    r3 = ad_us.gamma_ns_qcd(3, PID_NSM, c, nf3, var)
+    assert_allclose(r3, np.zeros_like(r3), atol=2e-4)
+
+    r3 = ad_us.gamma_ns_qcd(3, PID_NSV, c, nf3, var)
+    assert_allclose(r3, np.zeros_like(r3), atol=8e-4)
+
+    for v in range(3):
+        var_nsm = (0, v, 0)
+        r4 = ad_us.gamma_ns_qcd(4, PID_NSM, c, nf5, var_nsm)
+        assert_allclose(r4[3], nsm_refs[v], atol=6e-6)
+
+        var_nsv = (0, 0, v)
+        r4 = ad_us.gamma_ns_qcd(4, PID_NSV, c, nf5, var_nsv)
+        assert_allclose(r4[3], nsm_refs[v] + nss_refs[v], atol=1e-5)
+
+    r4 = ad_us.gamma_ns_qcd(4, PID_NSP, c, nf3, var)
+    assert np.any(np.abs(r4) ** 2 > 1e-12), "expected a non-trivial N3LO entry"
+
+
+def test_gamma_singlet_qcd():
+    nf = 5
+    c = Cache(2.0 + 0.0j)
+    quark_refs = [0.053441, 0.225674, -0.118792]
+    gluon_refs = [-0.0300842, 0.283004, -0.343172]
+
+    for imod in range(3):
+        var = (imod, imod, imod, imod)
+        g4 = ad_us.gamma_singlet_qcd(4, c, nf, var)
+
+        qq, qg = g4[3, 0, 0], g4[3, 0, 1]
+        gq, gg = g4[3, 1, 0], g4[3, 1, 1]
+
+        assert_allclose(qq + gq, quark_refs[imod], atol=2e-6)
+        assert_allclose(qg + gg, gluon_refs[imod], atol=2e-5)
+
+
+def test_gamma_ns_qed():
+    nf = 3
+    var = (0, 0, 0)
+    c = Cache(1.0 + 0.0j)
+    nsm_pids = [PID_NSM_U, PID_NSM_D]
+
+    for pid in nsm_pids:
+        r1 = ad_us.gamma_ns_qed(1, 1, pid, c, nf, var)
+        assert_allclose(r1, np.zeros_like(r1), atol=1e-5)
+
+    r1 = ad_us.gamma_ns_qed(1, 1, PID_NSP_U, c, nf, var)
+    assert_allclose(r1[0, 0], 0.0, atol=1e-15)
+    assert_allclose(r1[0, 1], 0.0, atol=1e-5)
+
+    r1 = ad_us.gamma_ns_qed(1, 1, PID_NSP_D, c, nf, var)
+    assert_allclose(r1[0, 0], 0.0, atol=1e-15)
+    assert_allclose(r1[0, 1], 0.0, atol=1e-5)
+
+    for pid in nsm_pids:
+        r2 = ad_us.gamma_ns_qed(1, 2, pid, c, nf, var)
+        assert_allclose(r2, np.zeros_like(r2), atol=1e-5)
+
+    for pid in nsm_pids:
+        r2 = ad_us.gamma_ns_qed(2, 1, pid, c, nf, var)
+        assert_allclose(r2, np.zeros_like(r2), atol=1e-5)
+
+    for pid in nsm_pids:
+        r3 = ad_us.gamma_ns_qed(3, 1, pid, c, nf, var)
+        assert_allclose(r3, np.zeros_like(r3), atol=1e-3)
+
+
+def test_gamma_valence_qed():
+    nf = 3
+    var = (0, 0, 0)
+    c = Cache(2.0 + 0.0j)
+
+    g = ad_us.gamma_valence_qed(3, 2, c, nf, var)
+
+    assert_allclose(g[0, 0], np.zeros((2, 2)), atol=1e-15)
+
+    assert_allclose(g[3, 0, 0, 0], 459.646893789751, atol=1e-5)
+    assert_allclose(g[3, 0, 0, 1], 0.0, atol=1e-5)
+    assert_allclose(g[3, 0, 1, 0], 0.0, atol=1e-5)
+    assert_allclose(g[3, 0, 1, 1], 437.60340375, atol=1e-5)
+
+
+def test_gamma_singlet_qed():
+    nf = 3
+    var = (0, 0, 0, 0, 0, 0, 0)
+    c = Cache(2.0 + 0.0j)
+
+    g = ad_us.gamma_singlet_qed(3, 2, c, nf, var)
+
+    assert_allclose(g[0, 0], np.zeros((4, 4)), atol=1e-15)
+
+    ref = np.array(
+        [
+            [3.857918949669738, 0.0, -290.42193908689745, 0.0],
+            [0.0, 0.0, 0.0, 0.0],
+            [-3.859554320251334, 0.0, 290.4252052962147, 0.0],
+            [0.0, 0.0, 0.0, 448.0752570151872],
+        ]
+    )
+    assert_allclose(g[3, 0].real, ref, atol=1e-5)
+    assert_allclose(g[3, 0].imag, np.zeros((4, 4)), atol=1e-5)
diff --git a/crates/ekore_py/tests/test_constants.py b/crates/ekore_py/tests/test_constants.py
new file mode 100644
index 000000000..d5bca9e10
--- /dev/null
+++ b/crates/ekore_py/tests/test_constants.py
@@ -0,0 +1,16 @@
+import ekore_py
+
+
+def test_max_orders():
+    assert ekore_py.constants.MAX_ORDER_QCD == 4
+    assert ekore_py.constants.MAX_ORDER_QED == 2
+
+
+def test_pids():
+    assert ekore_py.constants.PID_NSP == 10101
+    assert ekore_py.constants.PID_NSM == 10201
+    assert ekore_py.constants.PID_NSV == 10200
+    assert ekore_py.constants.PID_NSP_U == 10102
+    assert ekore_py.constants.PID_NSP_D == 10103
+    assert ekore_py.constants.PID_NSM_U == 10202
+    assert ekore_py.constants.PID_NSM_D == 10203
diff --git a/crates/ekore_py/tests/test_ome_us.py b/crates/ekore_py/tests/test_ome_us.py
new file mode 100644
index 000000000..58e2e9b86
--- /dev/null
+++ b/crates/ekore_py/tests/test_ome_us.py
@@ -0,0 +1,64 @@
+import numpy as np
+import pytest
+from ekore_py import Cache, ome_us
+from numpy.testing import assert_allclose
+
+
+@pytest.mark.parametrize(
+    "order, expected_shape", [(0, (0, 3, 3)), (1, (1, 3, 3)), (2, (2, 3, 3))]
+)
+def test_a_singlet_shape(order, expected_shape):
+    c = Cache(0.0 + 0.0j)
+    a = ome_us.A_singlet(order, c, nf=3, L=0.0)
+    assert a.shape == expected_shape
+    assert a.dtype == np.complex128
+
+
+@pytest.mark.parametrize(
+    "order, expected_shape", [(0, (0, 2, 2)), (1, (1, 2, 2)), (2, (2, 2, 2))]
+)
+def test_a_non_singlet_shape(order, expected_shape):
+    c = Cache(0.0 + 0.0j)
+    a = ome_us.A_non_singlet(order, c, nf=3, L=0.0)
+    assert a.shape == expected_shape
+    assert a.dtype == np.complex128
+
+
+def test_a_non_singlet():
+    nf = 5
+    L = 0.0
+    c = Cache(1.0 + 0.0j)
+
+    a = ome_us.A_non_singlet(2, c, nf, L)
+
+    # LO
+    assert_allclose(a[0], np.zeros((2, 2)), atol=1e-14)
+    # NNLO
+    assert_allclose(a[1], np.zeros((2, 2)), atol=1e-14)
+
+
+def test_a_singlet():
+    nf = 5
+    L = 100.0
+    c = Cache(2.0 + 0.0j)
+
+    a = ome_us.A_singlet(2, c, nf, L)
+
+    # LO
+    assert_allclose(a[0, 0, 2] + a[0, 1, 2] + a[0, 2, 2], 0.0, atol=1e-10)
+    assert_allclose(a[0, 0, 0] + a[0, 1, 0] + a[0, 2, 0], 0.0, atol=1e-12)
+
+    # NNLO
+    assert_allclose(a[1, 0, 0] + a[1, 1, 0] + a[1, 2, 0], 0.0, atol=2e-6)
+    assert_allclose(a[1, 0, 1] + a[1, 1, 1] + a[1, 2, 1], 0.0, atol=1e-11)
+
+
+def test_guards_invalid_order_raises():
+    nf = 4
+    L = 0.0
+    c = Cache(1.234 + 0.0j)
+
+    with pytest.raises(ValueError):
+        ome_us.A_singlet(3, c, nf, L)
+    with pytest.raises(ValueError):
+        ome_us.A_non_singlet(3, c, nf, L)
diff --git a/crates/release.json b/crates/release.json
index 0b51a21aa..9fe99bfa8 100644
--- a/crates/release.json
+++ b/crates/release.json
@@ -1 +1 @@
-["dekoder", "ekore", "eko", "ekore_capi"]
+["dekoder", "ekore", "eko", "ekore_capi", "ekore_py"]
diff --git a/pyproject.toml b/pyproject.toml
index 22924b260..da181eab2 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -156,6 +156,8 @@ clippy = "cargo clippy --no-deps"
 rbib = { "shell" = "python crates/make_bib.py > crates/ekore/src/bib.rs" }
 build-capi = "cargo cinstall --release -p ekore_capi --destdir=./crates/ekore_capi/dist --prefix=/ --libdir=/lib --includedir=/include"
 ctest = "./crates/ekore_capi/tests/run_tests.sh"
+build-pyapi = "maturin develop --release --manifest-path crates/ekore_py/Cargo.toml"
+ptest = { "shell" = "cd crates/ekore_py && pytest" }
 
 [tool.pytest.ini_options]
 testpaths = ['tests/', 'benchmarks/']