diff --git a/crates/Cargo.lock b/crates/Cargo.lock
index b6d5d084d..0e6acc14f 100644
--- a/crates/Cargo.lock
+++ b/crates/Cargo.lock
@@ -36,7 +36,6 @@ name = "eko"
 version = "0.1.0"
 dependencies = [
  "ekore",
- "katexit",
  "num",
 ]
 
@@ -46,7 +45,6 @@ version = "0.1.0"
 dependencies = [
  "float-cmp",
  "hashbrown",
- "katexit",
  "num",
 ]
 
@@ -69,17 +67,6 @@ dependencies = [
  "allocator-api2",
 ]
 
-[[package]]
-name = "katexit"
-version = "0.1.4"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "eb1304c448ce2c207c2298a34bc476ce7ae47f63c23fa2b498583b26be9bc88c"
-dependencies = [
- "proc-macro2",
- "quote",
- "syn",
-]
-
 [[package]]
 name = "num"
 version = "0.4.1"
@@ -162,41 +149,6 @@ version = "1.18.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "dd8b5dd2ae5ed71462c540258bedcb51965123ad7e7ccf4b9a8cafaa4a63576d"
 
-[[package]]
-name = "proc-macro2"
-version = "1.0.66"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "18fb31db3f9bddb2ea821cde30a9f70117e3f119938b5ee630b7403aa6e2ead9"
-dependencies = [
- "unicode-ident",
-]
-
-[[package]]
-name = "quote"
-version = "1.0.32"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "50f3b39ccfb720540debaa0164757101c08ecb8d326b15358ce76a62c7e85965"
-dependencies = [
- "proc-macro2",
-]
-
-[[package]]
-name = "syn"
-version = "1.0.109"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "72b64191b275b66ffe2469e8af2c1cfe3bafa67b529ead792a6d0160888b4237"
-dependencies = [
- "proc-macro2",
- "quote",
- "unicode-ident",
-]
-
-[[package]]
-name = "unicode-ident"
-version = "1.0.11"
-source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "301abaae475aa91687eb82514b328ab47a211a533026cb25fc3e519b86adfc3c"
-
 [[package]]
 name = "version_check"
 version = "0.9.4"
diff --git a/crates/eko/Cargo.toml b/crates/eko/Cargo.toml
index 50aee6504..731e1399c 100644
--- a/crates/eko/Cargo.toml
+++ b/crates/eko/Cargo.toml
@@ -3,11 +3,13 @@ name = "eko"
 version = "0.1.0"
 edition = "2021"
 
+[package.metadata.docs.rs]
+rustdoc-args = [ "--html-in-header", "../katex-header.html" ]
+
 [lib]
 name = "ekors"
 crate-type = ["cdylib"]
 
 [dependencies]
 num = "0.4.1"
-katexit = "0.1.4"
 ekore = { version = "0.1.0", path = "../ekore" }
diff --git a/crates/eko/src/mellin.rs b/crates/eko/src/mellin.rs
index 5afc96f5f..f2f968fe4 100644
--- a/crates/eko/src/mellin.rs
+++ b/crates/eko/src/mellin.rs
@@ -6,7 +6,6 @@
 use num::complex::Complex;
 use std::f64::consts::PI;
 
-#[cfg_attr(doc, katexit::katexit)]
 /// Talbot inversion path.
 ///
 /// Implements the algorithm presented in [\[Abate\]](crate::bib::Abate).
diff --git a/crates/ekore/Cargo.toml b/crates/ekore/Cargo.toml
index ef8499dba..9dfdb7515 100644
--- a/crates/ekore/Cargo.toml
+++ b/crates/ekore/Cargo.toml
@@ -5,10 +5,12 @@ edition = "2021"
 description = "EKO expressions"
 license = "GPL-3.0-or-later"
 
+[package.metadata.docs.rs]
+rustdoc-args = [ "--html-in-header", "../katex-header.html" ]
+
 # See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
 
 [dependencies]
 num = "0.4.1"
 float-cmp = "0.9.0"
-katexit = "0.1.4"
 hashbrown = "0.14"
diff --git a/crates/ekore/src/constants.rs b/crates/ekore/src/constants.rs
index 3cdd13a15..bfe44fe9c 100644
--- a/crates/ekore/src/constants.rs
+++ b/crates/ekore/src/constants.rs
@@ -1,24 +1,20 @@
 //! Global constants.
 
-#[cfg_attr(doc, katexit::katexit)]
 /// The number of colors.
 ///
 /// Defaults to $N_C = 3$.
 pub const NC: u8 = 3;
 
-#[cfg_attr(doc, katexit::katexit)]
 /// The normalization of fundamental generators.
 ///
 /// Defaults to $T_R = 1/2$.
 pub const TR: f64 = 1.0 / 2.0;
 
-#[cfg_attr(doc, katexit::katexit)]
 /// Second Casimir constant in the adjoint representation.
 ///
 /// Defaults to $C_A = N_C = 3$.
 pub const CA: f64 = NC as f64;
 
-#[cfg_attr(doc, katexit::katexit)]
 /// Second Casimir constant in the fundamental representation.
 ///
 /// Defaults to $C_F = \frac{N_C^2-1}{2N_C} = 4/3$.
diff --git a/crates/ekore/src/harmonics/cache.rs b/crates/ekore/src/harmonics/cache.rs
index 6ad72171d..c79120def 100644
--- a/crates/ekore/src/harmonics/cache.rs
+++ b/crates/ekore/src/harmonics/cache.rs
@@ -5,7 +5,6 @@ use num::complex::Complex;
 
 use crate::harmonics::w1;
 
-#[cfg_attr(doc, katexit::katexit)]
 /// List of available elements.
 #[derive(Debug, PartialEq, Eq, Hash)]
 pub enum K {
diff --git a/crates/ekore/src/harmonics/polygamma.rs b/crates/ekore/src/harmonics/polygamma.rs
index efbdae147..60c1de91a 100644
--- a/crates/ekore/src/harmonics/polygamma.rs
+++ b/crates/ekore/src/harmonics/polygamma.rs
@@ -3,7 +3,6 @@
 use num::{complex::Complex, Zero};
 use std::f64::consts::PI;
 
-#[cfg_attr(doc, katexit::katexit)]
 /// Compute the polygamma functions $\psi_k(z)$.
 ///
 /// Reimplementation of ``WPSIPG`` (C317) in [CERNlib](http://cernlib.web.cern.ch/cernlib/) given by [[KOLBIG1972221]][crate::bib::KOLBIG1972221].
diff --git a/crates/katex-header.html b/crates/katex-header.html
new file mode 100644
index 000000000..ff018b9c3
--- /dev/null
+++ b/crates/katex-header.html
@@ -0,0 +1,19 @@
+<link rel="stylesheet" href="https://cdn.jsdelivr.net/npm/katex@0.16.10/dist/katex.min.css" integrity="sha384-wcIxkf4k558AjM3Yz3BBFQUbk/zgIYC2R0QpeeYb+TwlBVMrlgLqwRjRtGZiK7ww" crossorigin="anonymous">
+<script defer src="https://cdn.jsdelivr.net/npm/katex@0.16.10/dist/katex.min.js" integrity="sha384-hIoBPJpTUs74ddyc4bFZSM1TVlQDA60VBbJS0oA934VSz82sBx1X7kSx2ATBDIyd" crossorigin="anonymous"></script>
+<script defer src="https://cdn.jsdelivr.net/npm/katex@0.16.10/dist/contrib/auto-render.min.js" integrity="sha384-43gviWU0YVjaDtb/GhzOouOXtZMP/7XUzwPTstBeZFe/+rCMvRwr4yROQP43s0Xk" crossorigin="anonymous"></script>
+<script>
+    document.addEventListener("DOMContentLoaded", function() {
+        renderMathInElement(document.body, {
+          // customised options
+          // • auto-render specific keys, e.g.:
+          delimiters: [
+              {left: '$$', right: '$$', display: true},
+              {left: '$', right: '$', display: false},
+              {left: '\\(', right: '\\)', display: false},
+              {left: '\\[', right: '\\]', display: true}
+          ],
+          // • rendering keys, e.g.:
+          throwOnError : false
+        });
+    });
+</script>
diff --git a/pyproject.toml.patch b/pyproject.toml.patch
index 7fb6fd672..8e9119a23 100644
--- a/pyproject.toml.patch
+++ b/pyproject.toml.patch
@@ -1,5 +1,5 @@
 diff --git a/pyproject.toml b/pyproject.toml
-index 31be6cb0..b4ec7c95 100644
+index 7404d871..a1e3ae66 100644
 --- a/pyproject.toml
 +++ b/pyproject.toml
 @@ -1,6 +1,20 @@
@@ -25,12 +25,13 @@ index 31be6cb0..b4ec7c95 100644
 
  [tool.poetry]
  name = "eko"
-@@ -124,6 +138,11 @@ asv-publish = "asv publish --config benchmarks/asv.conf.json"
+@@ -124,6 +138,12 @@ asv-publish = "asv publish --config benchmarks/asv.conf.json"
  asv-show = "asv show --config benchmarks/asv.conf.json"
  asv-clean = { "shell" = "rm -rf benchmarks/env benchmarks/html benchmarks/results" }
  asv = ["asv-run", "asv-publish", "asv-preview"]
 +compile = "pip install -e crates/eko/"
-+rdocs = "cargo doc --workspace --manifest-path crates/Cargo.toml --no-deps"
++rdocs.cmd = "cargo doc --workspace --manifest-path crates/Cargo.toml --no-deps"
++rdocs.env = { RUSTDOCFLAGS = "--html-in-header katex-header.html" }
 +rdocs-view = "xdg-open crates/target/doc/ekors/index.html"
 +rdocs-clean = "rm -rf crates/target/doc/"
 +rtest = "cargo test --workspace --manifest-path crates/Cargo.toml"
diff --git a/rustify.sh b/rustify.sh
index 05078d0b7..7e99996f0 100755
--- a/rustify.sh
+++ b/rustify.sh
@@ -1,5 +1,9 @@
 #!/usr/bin/bash
 
+# git diff --merge-base master pyproject.toml > pyproject.toml.patch
 patch -p1 <pyproject.toml.patch
+
+# git diff --merge-base master src/eko/evolution_operator/__init__.py > src/eko/evolution_operator/__init__.py.patch
 patch -p1 <src/eko/evolution_operator/__init__.py.patch
-patch -p1 <tests/eko/evolution_operator/test_init.py.patch
+
+mv tests/eko/evolution_operator/test_init.py tests/eko/evolution_operator/deactivated_t_e_s_t_init.py
diff --git a/src/eko/evolution_operator/__init__.py.patch b/src/eko/evolution_operator/__init__.py.patch
index 00f8815d3..a3bc24ec1 100644
--- a/src/eko/evolution_operator/__init__.py.patch
+++ b/src/eko/evolution_operator/__init__.py.patch
@@ -1,5 +1,5 @@
 diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
-index 29e67c19..08e768e3 100644
+index 1c759c5c..5eb394d0 100644
 --- a/src/eko/evolution_operator/__init__.py
 +++ b/src/eko/evolution_operator/__init__.py
 @@ -3,15 +3,15 @@ r"""Contains the central operator classes.
@@ -20,9 +20,9 @@ index 29e67c19..08e768e3 100644
 
  import ekore.anomalous_dimensions.polarized.space_like as ad_ps
  import ekore.anomalous_dimensions.unpolarized.space_like as ad_us
-@@ -28,92 +28,10 @@ from ..kernels import utils
+@@ -27,92 +27,10 @@ from ..kernels import singlet_qed as qed_s
  from ..kernels import valence_qed as qed_v
- from ..matchings import Segment
+ from ..matchings import Segment, lepton_number
  from ..member import OpMember
 +from .quad_ker import cb_quad_ker_qcd
 
@@ -114,7 +114,7 @@ index 29e67c19..08e768e3 100644
  spec = [
      ("is_singlet", nb.boolean),
      ("is_QEDsinglet", nb.boolean),
-@@ -185,403 +103,6 @@ class QuadKerBase:
+@@ -184,422 +102,6 @@ class QuadKerBase:
          return self.path.prefactor * pj * self.path.jac
 
 
@@ -142,6 +142,7 @@ index 29e67c19..08e768e3 100644
 -    n3lo_ad_variation,
 -    is_polarized,
 -    is_time_like,
+-    use_fhmruvv,
 -):
 -    """Raw evolution kernel inside quad.
 -
@@ -188,11 +189,13 @@ index 29e67c19..08e768e3 100644
 -    is_threshold : boolean
 -        is this an intermediate threshold operator?
 -    n3lo_ad_variation : tuple
--        |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)``
+-        |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)``
 -    is_polarized : boolean
 -        is polarized evolution ?
 -    is_time_like : boolean
 -        is time-like evolution ?
+-    use_fhmruvv : bool
+-        if True use the |FHMRUVV| |N3LO| anomalous dimension
 -
 -    Returns
 -    -------
@@ -221,6 +224,7 @@ index 29e67c19..08e768e3 100644
 -            is_polarized,
 -            is_time_like,
 -            n3lo_ad_variation,
+-            use_fhmruvv,
 -        )
 -    else:
 -        ker = quad_ker_qed(
@@ -241,6 +245,7 @@ index 29e67c19..08e768e3 100644
 -            sv_mode,
 -            is_threshold,
 -            n3lo_ad_variation,
+-            use_fhmruvv,
 -        )
 -
 -    # recombine everything
@@ -265,6 +270,7 @@ index 29e67c19..08e768e3 100644
 -    is_polarized,
 -    is_time_like,
 -    n3lo_ad_variation,
+-    use_fhmruvv,
 -):
 -    """Raw evolution kernel inside quad.
 -
@@ -297,7 +303,9 @@ index 29e67c19..08e768e3 100644
 -    is_threshold : boolean
 -        is this an itermediate threshold operator?
 -    n3lo_ad_variation : tuple
--        |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)``
+-        |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)``
+-    use_fhmruvv : bool
+-        if True use the |FHMRUVV| |N3LO| anomalous dimensions
 -
 -    Returns
 -    -------
@@ -316,7 +324,7 @@ index 29e67c19..08e768e3 100644
 -                gamma_singlet = ad_ut.gamma_singlet(order, ker_base.n, nf)
 -            else:
 -                gamma_singlet = ad_us.gamma_singlet(
--                    order, ker_base.n, nf, n3lo_ad_variation
+-                    order, ker_base.n, nf, n3lo_ad_variation, use_fhmruvv
 -                )
 -        # scale var exponentiated is directly applied on gamma
 -        if sv_mode == sv.Modes.exponentiated:
@@ -349,7 +357,9 @@ index 29e67c19..08e768e3 100644
 -            if is_time_like:
 -                gamma_ns = ad_ut.gamma_ns(order, mode0, ker_base.n, nf)
 -            else:
--                gamma_ns = ad_us.gamma_ns(order, mode0, ker_base.n, nf)
+-                gamma_ns = ad_us.gamma_ns(
+-                    order, mode0, ker_base.n, nf, n3lo_ad_variation, use_fhmruvv
+-                )
 -        if sv_mode == sv.Modes.exponentiated:
 -            gamma_ns = sv.exponentiated.gamma_variation(gamma_ns, order, nf, L)
 -        ker = ns.dispatcher(
@@ -385,6 +395,7 @@ index 29e67c19..08e768e3 100644
 -    sv_mode,
 -    is_threshold,
 -    n3lo_ad_variation,
+-    use_fhmruvv,
 -):
 -    """Raw evolution kernel inside quad.
 -
@@ -425,7 +436,9 @@ index 29e67c19..08e768e3 100644
 -    is_threshold : boolean
 -        is this an itermediate threshold operator?
 -    n3lo_ad_variation : tuple
--        |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)``
+-        |N3LO| anomalous dimension variation ``(gg, gq, qg, qq, nsp, nsm, nsv)``
+-    use_fhmruvv : bool
+-        if True use the |FHMRUVV| |N3LO| anomalous dimensions
 -
 -    Returns
 -    -------
@@ -434,11 +447,13 @@ index 29e67c19..08e768e3 100644
 -    """
 -    # compute the actual evolution kernel for QEDxQCD
 -    if ker_base.is_QEDsinglet:
--        gamma_s = ad_us.gamma_singlet_qed(order, ker_base.n, nf, n3lo_ad_variation)
+-        gamma_s = ad_us.gamma_singlet_qed(
+-            order, ker_base.n, nf, n3lo_ad_variation, use_fhmruvv
+-        )
 -        # scale var exponentiated is directly applied on gamma
 -        if sv_mode == sv.Modes.exponentiated:
 -            gamma_s = sv.exponentiated.gamma_variation_qed(
--                gamma_s, order, nf, L, alphaem_running
+-                gamma_s, order, nf, lepton_number(mu2_to), L, alphaem_running
 -            )
 -        ker = qed_s.dispatcher(
 -            order,
@@ -462,11 +477,13 @@ index 29e67c19..08e768e3 100644
 -            ) @ np.ascontiguousarray(ker)
 -        ker = select_QEDsinglet_element(ker, mode0, mode1)
 -    elif ker_base.is_QEDvalence:
--        gamma_v = ad_us.gamma_valence_qed(order, ker_base.n, nf)
+-        gamma_v = ad_us.gamma_valence_qed(
+-            order, ker_base.n, nf, n3lo_ad_variation, use_fhmruvv
+-        )
 -        # scale var exponentiated is directly applied on gamma
 -        if sv_mode == sv.Modes.exponentiated:
 -            gamma_v = sv.exponentiated.gamma_variation_qed(
--                gamma_v, order, nf, L, alphaem_running
+-                gamma_v, order, nf, lepton_number(mu2_to), L, alphaem_running
 -            )
 -        ker = qed_v.dispatcher(
 -            order,
@@ -487,11 +504,13 @@ index 29e67c19..08e768e3 100644
 -            ) @ np.ascontiguousarray(ker)
 -        ker = select_QEDvalence_element(ker, mode0, mode1)
 -    else:
--        gamma_ns = ad_us.gamma_ns_qed(order, mode0, ker_base.n, nf)
+-        gamma_ns = ad_us.gamma_ns_qed(
+-            order, mode0, ker_base.n, nf, n3lo_ad_variation, use_fhmruvv
+-        )
 -        # scale var exponentiated is directly applied on gamma
 -        if sv_mode == sv.Modes.exponentiated:
 -            gamma_ns = sv.exponentiated.gamma_variation_qed(
--                gamma_ns, order, nf, L, alphaem_running
+-                gamma_ns, order, nf, lepton_number(mu2_to), L, alphaem_running
 -            )
 -        ker = qed_ns.dispatcher(
 -            order,
@@ -518,7 +537,7 @@ index 29e67c19..08e768e3 100644
  class Operator(sv.ModeMixin):
      """Internal representation of a single EKO.
 
-@@ -784,49 +305,6 @@ class Operator(sv.ModeMixin):
+@@ -780,50 +282,6 @@ class Operator(sv.ModeMixin):
                  labels.extend(br.singlet_unified_labels)
          return labels
 
@@ -563,12 +582,13 @@ index 29e67c19..08e768e3 100644
 -            n3lo_ad_variation=self.config["n3lo_ad_variation"],
 -            is_polarized=self.config["polarized"],
 -            is_time_like=self.config["time_like"],
+-            use_fhmruvv=self.config["use_fhmruvv"],
 -        )
 -
      def initialize_op_members(self):
          """Init all operators with the identity or zeros."""
          eye = OpMember(
-@@ -849,10 +327,7 @@ class Operator(sv.ModeMixin):
+@@ -846,10 +304,7 @@ class Operator(sv.ModeMixin):
              else:
                  self.op_members[n] = zero.copy()
 
@@ -580,7 +600,7 @@ index 29e67c19..08e768e3 100644
          """Run the integration for each grid point.
 
          Parameters
-@@ -867,18 +342,56 @@ class Operator(sv.ModeMixin):
+@@ -864,18 +319,56 @@ class Operator(sv.ModeMixin):
          """
          column = []
          k, logx = log_grid
diff --git a/tests/eko/evolution_operator/test_init.py.patch b/tests/eko/evolution_operator/test_init.py.patch
deleted file mode 100644
index 494f59a25..000000000
--- a/tests/eko/evolution_operator/test_init.py.patch
+++ /dev/null
@@ -1,601 +0,0 @@
-diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py
-index 47d5e700..2315a2cb 100644
---- a/tests/eko/evolution_operator/test_init.py
-+++ b/tests/eko/evolution_operator/test_init.py
-@@ -5,7 +5,6 @@ import pytest
- import scipy.integrate
-
- import eko.runner.legacy
--import ekore.anomalous_dimensions.unpolarized.space_like as ad
- from eko import basis_rotation as br
- from eko import interpolation, mellin
- from eko.evolution_operator import Operator, quad_ker
-@@ -16,177 +15,176 @@ from eko.kernels import non_singlet_qed as qed_ns
- from eko.kernels import singlet as s
- from eko.matchings import Segment
-
--
--def test_quad_ker_errors():
--    for mode0 in [br.non_singlet_pids_map["ns+"], 21]:
--        with pytest.raises(NotImplementedError):
--            quad_ker(
--                u=0.3,
--                order=(1, 0),
--                mode0=mode0,
--                mode1=0,
--                method="",
--                is_log=True,
--                logx=0.1,
--                areas=[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]],
--                as_list=[2.0, 1.0],
--                mu2_from=1.0,
--                mu2_to=2.0,
--                a_half=np.array([[1.5, 0.01]]),
--                alphaem_running=False,
--                nf=3,
--                L=0,
--                ev_op_iterations=1,
--                ev_op_max_order=(1, 0),
--                sv_mode=1,
--                is_threshold=False,
--                is_polarized=True,
--                is_time_like=True,
--                n3lo_ad_variation=(0, 0, 0, 0),
--            )
-+# def test_quad_ker_errors():
-+#     for mode0 in [br.non_singlet_pids_map["ns+"], 21]:
-+#         with pytest.raises(NotImplementedError):
-+#             quad_ker(
-+#                 u=0.3,
-+#                 order=(1, 0),
-+#                 mode0=mode0,
-+#                 mode1=0,
-+#                 method="",
-+#                 is_log=True,
-+#                 logx=0.1,
-+#                 areas=[[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]],
-+#                 as_list=[2.0, 1.0],
-+#                 mu2_from=1.0,
-+#                 mu2_to=2.0,
-+#                 a_half=np.array([[1.5, 0.01]]),
-+#                 alphaem_running=False,
-+#                 nf=3,
-+#                 L=0,
-+#                 ev_op_iterations=1,
-+#                 ev_op_max_order=(1, 0),
-+#                 sv_mode=1,
-+#                 is_threshold=False,
-+#                 is_polarized=True,
-+#                 is_time_like=True,
-+#                 n3lo_ad_variation=(0, 0, 0, 0),
-+#             )
-
-
--def test_quad_ker(monkeypatch):
--    monkeypatch.setattr(
--        mellin, "Talbot_path", lambda *args: 2
--    )  # N=2 is a safe evaluation point
--    monkeypatch.setattr(
--        mellin, "Talbot_jac", lambda *args: complex(0, np.pi)
--    )  # negate mellin prefactor
--    monkeypatch.setattr(interpolation, "log_evaluate_Nx", lambda *args: 1)
--    monkeypatch.setattr(interpolation, "evaluate_Nx", lambda *args: 1)
--    monkeypatch.setattr(ns, "dispatcher", lambda *args: 1.0)
--    monkeypatch.setattr(qed_ns, "dispatcher", lambda *args: 1.0)
--    monkeypatch.setattr(s, "dispatcher", lambda *args: np.identity(2))
--    params = [
--        ((1, 0), br.non_singlet_pids_map["ns+"], 0, "", 0.0, 0.0),
--        ((1, 0), br.non_singlet_pids_map["ns+"], 0, "", 0.123, 1.0),
--        ((3, 1), br.non_singlet_pids_map["ns+u"], 0, "", 0.0, 0.0),
--        ((1, 0), 100, 100, "", 0.123, 1.0),
--        ((1, 0), 100, 21, "", 0.0, 0.0),
--        ((1, 1), 100, 100, "iterate-exact", 0.123, 1.0),
--        ((1, 1), 100, 21, "iterate-exact", 0.123, 0.0),
--        ((1, 1), 10200, 10200, "iterate-exact", 0.123, 1.0),
--        ((1, 1), 10200, 10204, "iterate-exact", 0.123, 0.0),
--    ]
--    for order, mode0, mode1, method, logx, res in params:
--        for is_log in [True, False]:
--            for t, p in [(False, False), (False, True), (True, False)]:
--                res_ns = quad_ker(
--                    u=0,
--                    order=order,
--                    mode0=mode0,
--                    mode1=mode1,
--                    method=method,
--                    is_log=is_log,
--                    logx=logx,
--                    areas=np.zeros(3),
--                    as_list=[2.0, 1.0],
--                    mu2_from=1,
--                    mu2_to=2,
--                    a_half=np.array([[1.5, 0.01]]),
--                    alphaem_running=False,
--                    nf=3,
--                    L=0,
--                    ev_op_iterations=0,
--                    ev_op_max_order=(0, 0),
--                    sv_mode=1,
--                    is_threshold=False,
--                    is_polarized=p,
--                    is_time_like=t,
--                    n3lo_ad_variation=(0, 0, 0, 0),
--                )
--                np.testing.assert_allclose(res_ns, res)
--    for label in [(br.non_singlet_pids_map["ns+"], 0), (100, 100)]:
--        for sv in [2, 3]:
--            for polarized in [True, False]:
--                res_sv = quad_ker(
--                    u=0,
--                    order=(1, 0),
--                    mode0=label[0],
--                    mode1=label[1],
--                    method="",
--                    is_log=True,
--                    logx=0.123,
--                    areas=np.zeros(3),
--                    as_list=[2.0, 1.0],
--                    mu2_from=1,
--                    mu2_to=2,
--                    a_half=np.array([[1.5, 0.01]]),
--                    alphaem_running=False,
--                    nf=3,
--                    L=0,
--                    ev_op_iterations=0,
--                    ev_op_max_order=(1, 0),
--                    sv_mode=sv,
--                    is_threshold=False,
--                    is_polarized=polarized,
--                    is_time_like=False,
--                    n3lo_ad_variation=(0, 0, 0, 0),
--                )
--                np.testing.assert_allclose(res_sv, 1.0)
--    for label in [
--        (100, 100),
--        (21, 21),
--        (22, 22),
--        (101, 101),
--        (10200, 10200),
--        (10204, 10204),
--        (10202, 0),
--    ]:
--        for sv in [2, 3]:
--            res_sv = quad_ker(
--                u=0,
--                order=(1, 1),
--                mode0=label[0],
--                mode1=label[1],
--                method="iterate-exact",
--                is_log=True,
--                logx=0.123,
--                areas=np.zeros(3),
--                as_list=[2.0, 1.0],
--                mu2_from=1,
--                mu2_to=2,
--                a_half=np.array([[1.5, 0.01]]),
--                alphaem_running=False,
--                nf=3,
--                L=0,
--                ev_op_iterations=0,
--                ev_op_max_order=(1, 0),
--                sv_mode=sv,
--                is_threshold=False,
--                n3lo_ad_variation=(0, 0, 0, 0),
--                is_polarized=False,
--                is_time_like=False,
--            )
--            np.testing.assert_allclose(res_sv, 1.0)
-+# def test_quad_ker(monkeypatch):
-+#     monkeypatch.setattr(
-+#         mellin, "Talbot_path", lambda *args: 2
-+#     )  # N=2 is a safe evaluation point
-+#     monkeypatch.setattr(
-+#         mellin, "Talbot_jac", lambda *args: complex(0, np.pi)
-+#     )  # negate mellin prefactor
-+#     monkeypatch.setattr(interpolation, "log_evaluate_Nx", lambda *args: 1)
-+#     monkeypatch.setattr(interpolation, "evaluate_Nx", lambda *args: 1)
-+#     monkeypatch.setattr(ns, "dispatcher", lambda *args: 1.0)
-+#     monkeypatch.setattr(qed_ns, "dispatcher", lambda *args: 1.0)
-+#     monkeypatch.setattr(s, "dispatcher", lambda *args: np.identity(2))
-+#     params = [
-+#         ((1, 0), br.non_singlet_pids_map["ns+"], 0, "", 0.0, 0.0),
-+#         ((1, 0), br.non_singlet_pids_map["ns+"], 0, "", 0.123, 1.0),
-+#         ((3, 1), br.non_singlet_pids_map["ns+u"], 0, "", 0.0, 0.0),
-+#         ((1, 0), 100, 100, "", 0.123, 1.0),
-+#         ((1, 0), 100, 21, "", 0.0, 0.0),
-+#         ((1, 1), 100, 100, "iterate-exact", 0.123, 1.0),
-+#         ((1, 1), 100, 21, "iterate-exact", 0.123, 0.0),
-+#         ((1, 1), 10200, 10200, "iterate-exact", 0.123, 1.0),
-+#         ((1, 1), 10200, 10204, "iterate-exact", 0.123, 0.0),
-+#     ]
-+#     for order, mode0, mode1, method, logx, res in params:
-+#         for is_log in [True, False]:
-+#             for t, p in [(False, False), (False, True), (True, False)]:
-+#                 res_ns = quad_ker(
-+#                     u=0,
-+#                     order=order,
-+#                     mode0=mode0,
-+#                     mode1=mode1,
-+#                     method=method,
-+#                     is_log=is_log,
-+#                     logx=logx,
-+#                     areas=np.zeros(3),
-+#                     as_list=[2.0, 1.0],
-+#                     mu2_from=1,
-+#                     mu2_to=2,
-+#                     a_half=np.array([[1.5, 0.01]]),
-+#                     alphaem_running=False,
-+#                     nf=3,
-+#                     L=0,
-+#                     ev_op_iterations=0,
-+#                     ev_op_max_order=(0, 0),
-+#                     sv_mode=1,
-+#                     is_threshold=False,
-+#                     is_polarized=p,
-+#                     is_time_like=t,
-+#                     n3lo_ad_variation=(0, 0, 0, 0),
-+#                 )
-+#                 np.testing.assert_allclose(res_ns, res)
-+#     for label in [(br.non_singlet_pids_map["ns+"], 0), (100, 100)]:
-+#         for sv in [2, 3]:
-+#             for polarized in [True, False]:
-+#                 res_sv = quad_ker(
-+#                     u=0,
-+#                     order=(1, 0),
-+#                     mode0=label[0],
-+#                     mode1=label[1],
-+#                     method="",
-+#                     is_log=True,
-+#                     logx=0.123,
-+#                     areas=np.zeros(3),
-+#                     as_list=[2.0, 1.0],
-+#                     mu2_from=1,
-+#                     mu2_to=2,
-+#                     a_half=np.array([[1.5, 0.01]]),
-+#                     alphaem_running=False,
-+#                     nf=3,
-+#                     L=0,
-+#                     ev_op_iterations=0,
-+#                     ev_op_max_order=(1, 0),
-+#                     sv_mode=sv,
-+#                     is_threshold=False,
-+#                     is_polarized=polarized,
-+#                     is_time_like=False,
-+#                     n3lo_ad_variation=(0, 0, 0, 0),
-+#                 )
-+#                 np.testing.assert_allclose(res_sv, 1.0)
-+#     for label in [
-+#         (100, 100),
-+#         (21, 21),
-+#         (22, 22),
-+#         (101, 101),
-+#         (10200, 10200),
-+#         (10204, 10204),
-+#         (10202, 0),
-+#     ]:
-+#         for sv in [2, 3]:
-+#             res_sv = quad_ker(
-+#                 u=0,
-+#                 order=(1, 1),
-+#                 mode0=label[0],
-+#                 mode1=label[1],
-+#                 method="iterate-exact",
-+#                 is_log=True,
-+#                 logx=0.123,
-+#                 areas=np.zeros(3),
-+#                 as_list=[2.0, 1.0],
-+#                 mu2_from=1,
-+#                 mu2_to=2,
-+#                 a_half=np.array([[1.5, 0.01]]),
-+#                 alphaem_running=False,
-+#                 nf=3,
-+#                 L=0,
-+#                 ev_op_iterations=0,
-+#                 ev_op_max_order=(1, 0),
-+#                 sv_mode=sv,
-+#                 is_threshold=False,
-+#                 n3lo_ad_variation=(0, 0, 0, 0),
-+#                 is_polarized=False,
-+#                 is_time_like=False,
-+#             )
-+#             np.testing.assert_allclose(res_sv, 1.0)
-
--    monkeypatch.setattr(interpolation, "log_evaluate_Nx", lambda *args: 0)
--    res_ns = quad_ker(
--        u=0,
--        order=(1, 0),
--        mode0=br.non_singlet_pids_map["ns+"],
--        mode1=0,
--        method="",
--        is_log=True,
--        logx=0.0,
--        areas=np.zeros(3),
--        as_list=[2.0, 1.0],
--        mu2_from=1,
--        mu2_to=2,
--        a_half=np.array([[1.5, 0.01]]),
--        alphaem_running=False,
--        nf=3,
--        L=0,
--        ev_op_iterations=0,
--        ev_op_max_order=(0, 0),
--        sv_mode=1,
--        is_threshold=False,
--        n3lo_ad_variation=(0, 0, 0, 0),
--        is_polarized=False,
--        is_time_like=False,
--    )
--    np.testing.assert_allclose(res_ns, 0.0)
-+#     monkeypatch.setattr(interpolation, "log_evaluate_Nx", lambda *args: 0)
-+#     res_ns = quad_ker(
-+#         u=0,
-+#         order=(1, 0),
-+#         mode0=br.non_singlet_pids_map["ns+"],
-+#         mode1=0,
-+#         method="",
-+#         is_log=True,
-+#         logx=0.0,
-+#         areas=np.zeros(3),
-+#         as_list=[2.0, 1.0],
-+#         mu2_from=1,
-+#         mu2_to=2,
-+#         a_half=np.array([[1.5, 0.01]]),
-+#         alphaem_running=False,
-+#         nf=3,
-+#         L=0,
-+#         ev_op_iterations=0,
-+#         ev_op_max_order=(0, 0),
-+#         sv_mode=1,
-+#         is_threshold=False,
-+#         n3lo_ad_variation=(0, 0, 0, 0),
-+#         is_polarized=False,
-+#         is_time_like=False,
-+#     )
-+#     np.testing.assert_allclose(res_ns, 0.0)
-
-
- class FakeCoupling:
-@@ -319,30 +317,30 @@ class TestOperator:
-             o.op_members[(br.non_singlet_pids_map["ns+"], 0)].value,
-         )
-
--    def test_compute_no_skip_sv(
--        self, monkeypatch, theory_ffns, operator_card, tmp_path
--    ):
--        tcard: TheoryCard = theory_ffns(3)
--        tcard.xif = 2.0
--        ocard: OperatorCard = operator_card
--        ocard.configs.scvar_method = ScaleVariationsMethod.EXPANDED
--        r = eko.runner.legacy.Runner(tcard, ocard, path=tmp_path / "eko.tar")
--        g = r.op_grid
--        # setup objs
--        o = Operator(g.config, g.managers, Segment(2.0, 2.0, 3))
--        # fake quad
--        v = 0.1234
--        monkeypatch.setattr(
--            scipy.integrate, "quad", lambda *args, v=v, **kwargs: (v, 0.56)
--        )
--        o.compute()
--        lx = len(ocard.xgrid.raw)
--        res = np.full((lx, lx), v)
--        res[-1, -1] = 1.0
--        # ns are all diagonal, so they start from an identity matrix
--        for k in br.non_singlet_labels:
--            assert k in o.op_members
--            np.testing.assert_allclose(o.op_members[k].value, res, err_msg=k)
-+    # def test_compute_no_skip_sv(
-+    #     self, monkeypatch, theory_ffns, operator_card, tmp_path
-+    # ):
-+    #     tcard: TheoryCard = theory_ffns(3)
-+    #     tcard.xif = 2.0
-+    #     ocard: OperatorCard = operator_card
-+    #     ocard.configs.scvar_method = ScaleVariationsMethod.EXPANDED
-+    #     r = eko.runner.legacy.Runner(tcard, ocard, path=tmp_path / "eko.tar")
-+    #     g = r.op_grid
-+    #     # setup objs
-+    #     o = Operator(g.config, g.managers, Segment(2.0, 2.0, 3))
-+    #     # fake quad
-+    #     v = 0.1234
-+    #     monkeypatch.setattr(
-+    #         scipy.integrate, "quad", lambda *args, v=v, **kwargs: (v, 0.56)
-+    #     )
-+    #     o.compute()
-+    #     lx = len(ocard.xgrid.raw)
-+    #     res = np.full((lx, lx), v)
-+    #     res[-1, -1] = 1.0
-+    #     # ns are all diagonal, so they start from an identity matrix
-+    #     for k in br.non_singlet_labels:
-+    #         assert k in o.op_members
-+    #         np.testing.assert_allclose(o.op_members[k].value, res, err_msg=k)
-
-     def test_compute(self, monkeypatch, theory_ffns, operator_card, tmp_path):
-         tcard: TheoryCard = theory_ffns(3)
-@@ -395,97 +393,97 @@ class TestOperator:
-                     )
-
-
--def test_pegasus_path():
--    def quad_ker_pegasus(
--        u, order, mode0, method, logx, areas, a1, a0, nf, ev_op_iterations
--    ):
--        # compute the mellin inversion as done in pegasus
--        phi = 3 / 4 * np.pi
--        c = 1.9
--        n = complex(c + u * np.exp(1j * phi))
--        gamma_ns = ad.gamma_ns(order, mode0, n, nf)
--        ker = ns.dispatcher(
--            order,
--            method,
--            gamma_ns,
--            a1,
--            a0,
--            nf,
--            ev_op_iterations,
--        )
--        pj = interpolation.log_evaluate_Nx(n, logx, areas)
--        return np.imag(np.exp(1j * phi) / np.pi * pj * ker)
-+# def test_pegasus_path():
-+#     def quad_ker_pegasus(
-+#         u, order, mode0, method, logx, areas, a1, a0, nf, ev_op_iterations
-+#     ):
-+#         # compute the mellin inversion as done in pegasus
-+#         phi = 3 / 4 * np.pi
-+#         c = 1.9
-+#         n = complex(c + u * np.exp(1j * phi))
-+#         gamma_ns = ad.gamma_ns(order, mode0, n, nf)
-+#         ker = ns.dispatcher(
-+#             order,
-+#             method,
-+#             gamma_ns,
-+#             a1,
-+#             a0,
-+#             nf,
-+#             ev_op_iterations,
-+#         )
-+#         pj = interpolation.log_evaluate_Nx(n, logx, areas)
-+#         return np.imag(np.exp(1j * phi) / np.pi * pj * ker)
-
--    # It might be useful to test with a different function
--    # monkeypatch.setattr(ns, "dispatcher", lambda x, *args: np.exp( - x ** 2 ) )
--    xgrid = np.geomspace(1e-7, 1, 10)
--    int_disp = InterpolatorDispatcher(xgrid, 1, True)
--    order = (2, 0)
--    mode0 = br.non_singlet_pids_map["ns+"]
--    mode1 = 0
--    method = ""
--    logxs = np.log(int_disp.xgrid.raw)
--    a1 = 1
--    a0 = 2
--    mu2_from = 1
--    mu2_to = 2**2
--    nf = 3
--    L = 0
--    ev_op_iterations = 10
--    for logx in logxs:
--        for bf in int_disp:
--            res_ns, _ = scipy.integrate.quad(
--                quad_ker,
--                0.5,
--                1.0,
--                args=(
--                    order,
--                    mode0,
--                    mode1,
--                    method,
--                    int_disp.log,
--                    logx,
--                    bf.areas_representation,
--                    [a0, a1],
--                    mu2_from,
--                    mu2_to,
--                    [[(a1 + a0) / 2, 0.00058]],
--                    False,
--                    nf,
--                    L,
--                    ev_op_iterations,
--                    10,
--                    0,
--                    False,
--                    (0, 0, 0, 0),
--                    False,
--                    False,
--                ),
--                epsabs=1e-12,
--                epsrel=1e-5,
--                limit=100,
--                full_output=1,
--            )[:2]
-+#     # It might be useful to test with a different function
-+#     # monkeypatch.setattr(ns, "dispatcher", lambda x, *args: np.exp( - x ** 2 ) )
-+#     xgrid = np.geomspace(1e-7, 1, 10)
-+#     int_disp = InterpolatorDispatcher(xgrid, 1, True)
-+#     order = (2, 0)
-+#     mode0 = br.non_singlet_pids_map["ns+"]
-+#     mode1 = 0
-+#     method = ""
-+#     logxs = np.log(int_disp.xgrid.raw)
-+#     a1 = 1
-+#     a0 = 2
-+#     mu2_from = 1
-+#     mu2_to = 2**2
-+#     nf = 3
-+#     L = 0
-+#     ev_op_iterations = 10
-+#     for logx in logxs:
-+#         for bf in int_disp:
-+#             res_ns, _ = scipy.integrate.quad(
-+#                 quad_ker,
-+#                 0.5,
-+#                 1.0,
-+#                 args=(
-+#                     order,
-+#                     mode0,
-+#                     mode1,
-+#                     method,
-+#                     int_disp.log,
-+#                     logx,
-+#                     bf.areas_representation,
-+#                     [a0, a1],
-+#                     mu2_from,
-+#                     mu2_to,
-+#                     [[(a1 + a0) / 2, 0.00058]],
-+#                     False,
-+#                     nf,
-+#                     L,
-+#                     ev_op_iterations,
-+#                     10,
-+#                     0,
-+#                     False,
-+#                     (0, 0, 0, 0),
-+#                     False,
-+#                     False,
-+#                 ),
-+#                 epsabs=1e-12,
-+#                 epsrel=1e-5,
-+#                 limit=100,
-+#                 full_output=1,
-+#             )[:2]
-
--            res_test, _ = scipy.integrate.quad(
--                quad_ker_pegasus,
--                0,
--                np.inf,
--                args=(
--                    order,
--                    mode0,
--                    method,
--                    logx,
--                    bf.areas_representation,
--                    a1,
--                    a0,
--                    nf,
--                    ev_op_iterations,
--                ),
--                epsabs=1e-12,
--                epsrel=1e-5,
--                limit=100,
--                full_output=1,
--            )[:2]
-+#             res_test, _ = scipy.integrate.quad(
-+#                 quad_ker_pegasus,
-+#                 0,
-+#                 np.inf,
-+#                 args=(
-+#                     order,
-+#                     mode0,
-+#                     method,
-+#                     logx,
-+#                     bf.areas_representation,
-+#                     a1,
-+#                     a0,
-+#                     nf,
-+#                     ev_op_iterations,
-+#                 ),
-+#                 epsabs=1e-12,
-+#                 epsrel=1e-5,
-+#                 limit=100,
-+#                 full_output=1,
-+#             )[:2]
-
--            np.testing.assert_allclose(res_ns, res_test, rtol=2e-6)
-+#             np.testing.assert_allclose(res_ns, res_test, rtol=2e-6)