diff --git a/crates/eko/src/lib.rs b/crates/eko/src/lib.rs
index 152db2d56..958cfc3fb 100644
--- a/crates/eko/src/lib.rs
+++ b/crates/eko/src/lib.rs
@@ -20,106 +20,68 @@ pub unsafe extern "C" fn rust_quad_ker_qcd(u: f64, rargs: *mut c_void) -> f64 {
     let mut c = Cache::new(path.n());
     let mut re = Vec::<f64>::new();
     let mut im = Vec::<f64>::new();
-    if is_singlet {
-        let res = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qcd(
-            args.order_qcd,
-            &mut c,
-            args.nf,
-        );
-        for gamma_s in res.iter().take(args.order_qcd) {
-            for col in gamma_s.iter().take(2) {
-                for el in col.iter().take(2) {
-                    re.push(el.re);
-                    im.push(el.im);
+
+    if args.is_ome {
+        if is_singlet {
+            let res = ekore::operator_matrix_elements::unpolarized::spacelike::A_singlet(
+                args.order_qcd,
+                &mut c,
+                args.nf,
+                args.L,
+            );
+            for aS in res.iter().take(args.order_qcd) {
+                for col in aS.iter().take(3) {
+                    for el in col.iter().take(3) {
+                        re.push(el.re);
+                        im.push(el.im);
+                    }
                 }
             }
-        }
-    } else {
-        let res = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_ns_qcd(
-            args.order_qcd,
-            args.mode0,
-            &mut c,
-            args.nf,
-        );
-        for el in res.iter().take(args.order_qcd) {
-            re.push(el.re);
-            im.push(el.im);
-        }
-    }
-    // pass on
-    (args.py)(
-        re.as_ptr(),
-        im.as_ptr(),
-        c.n.re,
-        c.n.im,
-        jac.re,
-        jac.im,
-        args.order_qcd,
-        is_singlet,
-        args.mode0,
-        args.mode1,
-        args.nf,
-        args.is_log,
-        args.logx,
-        args.areas,
-        args.areas_x,
-        args.areas_y,
-        args.L,
-        args.method_num,
-        args.as1,
-        args.as0,
-        args.ev_op_iterations,
-        args.ev_op_max_order_qcd,
-        args.sv_mode_num,
-        args.is_threshold,
-    )
-}
-
-/// QCD intergration kernel for OME inside quad.
-///
-/// # Safety
-/// This is the connection from Python, so we don't know what is on the other side.
-#[no_mangle]
-pub unsafe extern "C" fn rust_quad_ker_ome(u: f64, rargs: *mut c_void) -> f64 {
-    let args = *(rargs as *mut QuadQCDargs);
-    let is_singlet = (100 == args.mode0) || (21 == args.mode0) || (90 == args.mode0);
-    // prepare gamma
-    let path = mellin::TalbotPath::new(u, args.logx, is_singlet);
-    let jac = path.jac() * path.prefactor();
-    let mut c = Cache::new(path.n());
-    let mut re = Vec::<f64>::new();
-    let mut im = Vec::<f64>::new();
-    if is_singlet {
-        let res = ekore::operator_matrix_elements::unpolarized::spacelike::A_singlet(
-            args.order_qcd,
-            &mut c,
-            args.nf,
-            args.L,
-        );
-        for aS in res.iter().take(args.order_qcd) {
-            for col in aS.iter().take(3) {
-                for el in col.iter().take(3) {
-                    re.push(el.re);
-                    im.push(el.im);
+        } else {
+            let res = ekore::operator_matrix_elements::unpolarized::spacelike::A_non_singlet(
+                args.order_qcd,
+                &mut c,
+                args.nf,
+                args.L,
+            );
+            for anS in res.iter().take(args.order_qcd) {
+                for col in anS.iter().take(2) {
+                    for el in col.iter().take(2) {
+                        re.push(el.re);
+                        im.push(el.im);
+                    }
                 }
             }
         }
     } else {
-        let res = ekore::operator_matrix_elements::unpolarized::spacelike::A_non_singlet(
-            args.order_qcd,
-            &mut c,
-            args.nf,
-            args.L,
-        );
-        for anS in res.iter().take(args.order_qcd) {
-            for col in anS.iter().take(2) {
-                for el in col.iter().take(2) {
-                    re.push(el.re);
-                    im.push(el.im);
+        if is_singlet {
+            let res = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_singlet_qcd(
+                args.order_qcd,
+                &mut c,
+                args.nf,
+            );
+            for gamma_s in res.iter().take(args.order_qcd) {
+                for col in gamma_s.iter().take(2) {
+                    for el in col.iter().take(2) {
+                        re.push(el.re);
+                        im.push(el.im);
+                    }
                 }
             }
+        } else {
+            let res = ekore::anomalous_dimensions::unpolarized::spacelike::gamma_ns_qcd(
+                args.order_qcd,
+                args.mode0,
+                &mut c,
+                args.nf,
+            );
+            for el in res.iter().take(args.order_qcd) {
+                re.push(el.re);
+                im.push(el.im);
+            }
         }
     }
+
     // pass on
     (args.py)(
         re.as_ptr(),
@@ -146,6 +108,7 @@ pub unsafe extern "C" fn rust_quad_ker_ome(u: f64, rargs: *mut c_void) -> f64 {
         args.ev_op_max_order_qcd,
         args.sv_mode_num,
         args.is_threshold,
+        args.Lsv,
     )
 }
 
@@ -175,6 +138,7 @@ type PyQuadKerQCDT = unsafe extern "C" fn(
     u8,
     u8,
     bool,
+    f64,
 ) -> f64;
 
 /// Additional integration parameters
@@ -202,6 +166,8 @@ pub struct QuadQCDargs {
     pub ev_op_max_order_qcd: u8,
     pub sv_mode_num: u8,
     pub is_threshold: bool,
+    pub is_ome: bool,
+    pub Lsv: f64,
 }
 
 /// Empty placeholder function for python callback.
@@ -233,6 +199,7 @@ pub unsafe extern "C" fn my_py(
     _ev_op_max_order_qcd: u8,
     _sv_mode_num: u8,
     _is_threshold: bool,
+    _lsv: f64,
 ) -> f64 {
     0.
 }
@@ -267,5 +234,7 @@ pub unsafe extern "C" fn empty_qcd_args() -> QuadQCDargs {
         ev_op_max_order_qcd: 0,
         sv_mode_num: 0,
         is_threshold: false,
+        is_ome: false,
+        Lsv: 0.,
     }
 }
diff --git a/src/eko/evolution_operator/operator_matrix_element.py b/src/eko/evolution_operator/operator_matrix_element.py
index 202fd18d9..f76db16a5 100644
--- a/src/eko/evolution_operator/operator_matrix_element.py
+++ b/src/eko/evolution_operator/operator_matrix_element.py
@@ -1,9 +1,9 @@
 """The |OME| for the non-trivial matching conditions in the |VFNS| evolution."""
 
 import copy
-import functools
 import logging
 
+import ekors
 import numba as nb
 import numpy as np
 
@@ -15,7 +15,8 @@
 from .. import scale_variations as sv
 from ..io.types import InversionMethod
 from ..matchings import Segment
-from . import Operator, QuadKerBase
+from . import Operator
+from .quad_ker import cb_quad_ker_ome
 
 logger = logging.getLogger(__name__)
 
@@ -73,106 +74,6 @@ def build_ome(A, matching_order, a_s, backward_method):
     return ome
 
 
-@nb.njit(cache=True)
-def quad_ker(
-    u,
-    order,
-    mode0,
-    mode1,
-    is_log,
-    logx,
-    areas,
-    a_s,
-    nf,
-    L,
-    sv_mode,
-    Lsv,
-    backward_method,
-    is_msbar,
-    is_polarized,
-    is_time_like,
-):
-    r"""Raw kernel inside quad.
-
-    Parameters
-    ----------
-    u : float
-        quad argument
-    order : tuple(int,int)
-        perturbation matching order
-    mode0 : int
-        pid for first element in the singlet sector
-    mode1 : int
-        pid for second element in the singlet sector
-    is_log : boolean
-        logarithmic interpolation
-    logx : float
-        Mellin inversion point
-    areas : tuple
-        basis function configuration
-    a_s : float
-        strong coupling, needed only for the exact inverse
-    nf: int
-        number of active flavor below threshold
-    L : float
-        :math:``\ln(\mu_F^2 / m_h^2)``
-    backward_method : InversionMethod or None
-        empty or method for inverting the matching condition (exact or expanded)
-    is_msbar: bool
-        add the |MSbar| contribution
-    is_polarized : boolean
-        is polarized evolution ?
-    is_time_like : boolean
-        is time-like evolution ?
-
-    Returns
-    -------
-    ker : float
-        evaluated integration kernel
-
-    """
-    ker_base = QuadKerBase(u, is_log, logx, mode0)
-    integrand = ker_base.integrand(areas)
-    if integrand == 0.0:
-        return 0.0
-    # compute the ome
-    if ker_base.is_singlet or ker_base.is_QEDsinglet:
-        indices = {21: 0, 100: 1, 90: 2}
-        if is_polarized:
-            if is_time_like:
-                raise NotImplementedError("Polarized, time-like is not implemented")
-            A = ome_ps.A_singlet(order, ker_base.n, nf, L)
-        else:
-            if is_time_like:
-                A = ome_ut.A_singlet(order, ker_base.n, L)
-            else:
-                A = ome_us.A_singlet(order, ker_base.n, nf, L, is_msbar)
-    else:
-        indices = {200: 0, 91: 1}
-        if is_polarized:
-            if is_time_like:
-                raise NotImplementedError("Polarized, time-like is not implemented")
-            A = ome_ps.A_non_singlet(order, ker_base.n, L)
-        else:
-            if is_time_like:
-                A = ome_ut.A_non_singlet(order, ker_base.n, L)
-            else:
-                A = ome_us.A_non_singlet(order, ker_base.n, nf, L)
-
-    # correct for scale variations
-    if sv_mode == sv.Modes.exponentiated:
-        A = sv.exponentiated.gamma_variation(A, order, nf, Lsv)
-
-    # build the expansion in alpha_s depending on the strategy
-    ker = build_ome(A, order, a_s, backward_method)
-
-    # select the needed matrix element
-    ker = ker[indices[mode0], indices[mode1]]
-
-    # recombine everything
-    return np.real(ker * integrand)
-
-
 class OperatorMatrixElement(Operator):
     r"""
     Internal representation of a single |OME|.
@@ -267,41 +168,89 @@ def labels(self):
                 )
         return labels
 
-    def quad_ker(self, label, logx, areas):
-        """Return partially initialized integrand function.
+    def run_op_integration(self, log_grid):
+        """Run the integration for each grid point.
 
         Parameters
         ----------
-        label: tuple
-            operator element pids
-        logx: float
-            Mellin inversion point
-        areas : tuple
-            basis function configuration
+        log_grid : tuple(k, logx)
+            log grid point with relative index
 
         Returns
         -------
-        functools.partial
-            partially initialized integration kernel
+        list
+            computed operators at the give grid point
         """
-        return functools.partial(
-            quad_ker,
-            order=self.order,
-            mode0=label[0],
-            mode1=label[1],
-            is_log=self.int_disp.log,
-            logx=logx,
-            areas=areas,
-            a_s=self.a_s,
-            nf=self.nf,
-            L=self.L,
-            sv_mode=self.sv_mode,
-            Lsv=np.log(self.xif2),
-            backward_method=self.backward_method,
-            is_msbar=self.is_msbar,
-            is_polarized=self.config["polarized"],
-            is_time_like=self.config["time_like"],
+        column = []
+        k, logx = log_grid
+        # call(!) self.labels only once
+        labels = self.labels
+        start_time = time.perf_counter()
+        # start preparing C arguments
+        cfg = ekors.lib.empty_qcd_args()
+        cfg.order_qcd = self.order[0]
+        cfg.is_polarized = self.config["polarized"]
+        cfg.is_time_like = self.config["time_like"]
+        cfg.nf = self.nf
+        cfg.py = ekors.ffi.cast("void *", cb_quad_ker_ome.address)
+        cfg.is_log = self.int_disp.log
+        cfg.logx = logx
+        cfg.L = self.L
+        # cfg.method_num = 1
+        cfg.as1 = self.as_list[1]
+        cfg.as0 = self.as_list[0]
+        # cfg.ev_op_iterations = self.config["ev_op_iterations"]
+        # cfg.ev_op_max_order_qcd = self.config["ev_op_max_order"][0]
+        # cfg.sv_mode_num = 1
+        # cfg.is_threshold = self.is_threshold
+        cfg.Lsv = np.log(self.xif2)
+
+        # iterate basis functions
+        for l, bf in enumerate(self.int_disp):
+            if k == l and l == self.grid_size - 1:
+                continue
+            # add emtpy labels with 0s
+            if bf.is_below_x(np.exp(logx)):
+                column.append({label: (0.0, 0.0) for label in labels})
+                continue
+            temp_dict = {}
+            # prepare areas for C
+            curareas = bf.areas_representation
+            areas_len = curareas.shape[0] * curareas.shape[1]
+            # force the variable in scope
+            areas_ffi = ekors.ffi.new(
+                f"double[{areas_len}]", curareas.flatten().tolist()
+            )
+            cfg.areas = areas_ffi
+            cfg.areas_x = curareas.shape[0]
+            cfg.areas_y = curareas.shape[1]
+            # iterate sectors
+            for label in labels:
+                cfg.mode0 = label[0]
+                cfg.mode1 = label[1]
+                # construct the low level object
+                func = LowLevelCallable(
+                    ekors.lib.rust_quad_ker_ome, ekors.ffi.addressof(cfg)
+                )
+                res = integrate.quad(
+                    func,
+                    0.5,
+                    1.0 - self._mellin_cut,
+                    epsabs=1e-12,
+                    epsrel=1e-5,
+                    limit=100,
+                    full_output=1,
+                )
+                temp_dict[label] = res[:2]
+            column.append(temp_dict)
+        logger.info(
+            "%s: computing operators - %u/%u took: %6f s",
+            self.log_label,
+            k + 1,
+            self.grid_size,
+            (time.perf_counter() - start_time),
         )
+        return column
 
     @property
     def a_s(self):
diff --git a/src/eko/evolution_operator/operator_matrix_element_rust.py b/src/eko/evolution_operator/operator_matrix_element_rust.py
deleted file mode 100644
index 4cf1bd20d..000000000
--- a/src/eko/evolution_operator/operator_matrix_element_rust.py
+++ /dev/null
@@ -1,421 +0,0 @@
-"""The |OME| for the non-trivial matching conditions in the |VFNS| evolution."""
-
-import copy
-import logging
-
-import ekors
-import numba as nb
-import numpy as np
-
-import ekore.operator_matrix_elements.polarized.space_like as ome_ps
-import ekore.operator_matrix_elements.unpolarized.space_like as ome_us
-import ekore.operator_matrix_elements.unpolarized.time_like as ome_ut
-
-from .. import basis_rotation as br
-from .. import scale_variations as sv
-from ..io.types import InversionMethod
-from ..matchings import Segment
-from . import Operator
-from .quad_ker_ome import cb_quad_ker_ome
-
-logger = logging.getLogger(__name__)
-
-
-# @nb.njit(cache=True)
-# def build_ome(A, matching_order, a_s, backward_method):
-#     r"""Construct the matching expansion in :math:`a_s` with the appropriate method.
-
-#     Parameters
-#     ----------
-#     A : numpy.ndarray
-#         list of |OME|
-#     matching_order : tuple(int,int)
-#         perturbation matching order
-#     a_s : float
-#         strong coupling, needed only for the exact inverse
-#     backward_method : InversionMethod or None
-#         empty or method for inverting the matching condition (exact or expanded)
-
-#     Returns
-#     -------
-#     ome : numpy.ndarray
-#         matching operator matrix
-
-#     """
-#     # to get the inverse one can use this FORM snippet
-#     # Symbol a;
-#     # NTensor c,d,e;
-#     # Local x=-(a*c+a**2* d + a**3 * e);
-#     # Local bi = 1+x+x**2+x**3;
-#     # Print;
-#     # .end
-#     ome = np.eye(len(A[0]), dtype=np.complex_)
-#     A = A[:, :, :]
-#     A = np.ascontiguousarray(A)
-#     if backward_method is InversionMethod.EXPANDED:
-#         # expended inverse
-#         if matching_order[0] >= 1:
-#             ome -= a_s * A[0]
-#         if matching_order[0] >= 2:
-#             ome += a_s**2 * (-A[1] + A[0] @ A[0])
-#         if matching_order[0] >= 3:
-#             ome += a_s**3 * (-A[2] + A[0] @ A[1] + A[1] @ A[0] - A[0] @ A[0] @ A[0])
-#     else:
-#         # forward or exact inverse
-#         if matching_order[0] >= 1:
-#             ome += a_s * A[0]
-#         if matching_order[0] >= 2:
-#             ome += a_s**2 * A[1]
-#         if matching_order[0] >= 3:
-#             ome += a_s**3 * A[2]
-#         # need inverse exact ?  so add the missing pieces
-#         if backward_method is InversionMethod.EXACT:
-#             ome = np.linalg.inv(ome)
-#     return ome
-
-# @nb.njit(cache=True)
-# def quad_ker(
-#     u,
-#     order,
-#     mode0,
-#     mode1,
-#     is_log,
-#     logx,
-#     areas,
-#     a_s,
-#     nf,
-#     L,
-#     sv_mode,
-#     Lsv,
-#     backward_method,
-#     is_msbar,
-#     is_polarized,
-#     is_time_like,
-# ):
-#     r"""Raw kernel inside quad.
-
-#     Parameters
-#     ----------
-#     u : float
-#         quad argument
-#     order : tuple(int,int)
-#         perturbation matching order
-#     mode0 : int
-#         pid for first element in the singlet sector
-#     mode1 : int
-#         pid for second element in the singlet sector
-#     is_log : boolean
-#         logarithmic interpolation
-#     logx : float
-#         Mellin inversion point
-#     areas : tuple
-#         basis function configuration
-#     a_s : float
-#         strong coupling, needed only for the exact inverse
-#     nf: int
-#         number of active flavor below threshold
-#     L : float
-#         :math:``\ln(\mu_F^2 / m_h^2)``
-#     backward_method : InversionMethod or None
-#         empty or method for inverting the matching condition (exact or expanded)
-#     is_msbar: bool
-#         add the |MSbar| contribution
-#     is_polarized : boolean
-#         is polarized evolution ?
-#     is_time_like : boolean
-#         is time-like evolution ?
-
-#     Returns
-#     -------
-#     ker : float
-#         evaluated integration kernel
-
-#     """
-#     ker_base = QuadKerBase(u, is_log, logx, mode0)
-#     integrand = ker_base.integrand(areas)
-#     if integrand == 0.0:
-#         return 0.0
-#     # compute the ome
-#     if ker_base.is_singlet or ker_base.is_QEDsinglet:
-#         indices = {21: 0, 100: 1, 90: 2}
-#         if is_polarized:
-#             if is_time_like:
-#                 raise NotImplementedError("Polarized, time-like is not implemented")
-#             A = ome_ps.A_singlet(order, ker_base.n, nf, L)
-#         else:
-#             if is_time_like:
-#                 A = ome_ut.A_singlet(order, ker_base.n, L)
-#             else:
-#                 A = ome_us.A_singlet(order, ker_base.n, nf, L, is_msbar)
-#     else:
-#         indices = {200: 0, 91: 1}
-#         if is_polarized:
-#             if is_time_like:
-#                 raise NotImplementedError("Polarized, time-like is not implemented")
-#             A = ome_ps.A_non_singlet(order, ker_base.n, L)
-#         else:
-#             if is_time_like:
-#                 A = ome_ut.A_non_singlet(order, ker_base.n, L)
-#             else:
-#                 A = ome_us.A_non_singlet(order, ker_base.n, nf, L)
-
-#     # correct for scale variations
-#     if sv_mode == sv.Modes.exponentiated:
-#         A = sv.exponentiated.gamma_variation(A, order, nf, Lsv)
-
-#     # build the expansion in alpha_s depending on the strategy
-#     ker = build_ome(A, order, a_s, backward_method)
-
-#     # select the needed matrix element
-#     ker = ker[indices[mode0], indices[mode1]]
-
-#     # recombine everything
-#     return np.real(ker * integrand)
-
-
-# This has to be modified by implementing a run_op_integration function like done in __init__.py/Operator
-
-
-class OperatorMatrixElement(Operator):
-    r"""
-    Internal representation of a single |OME|.
-
-    The actual matrices are computed upon calling :meth:`compute`.
-
-    Parameters
-    ----------
-    config : dict
-        configuration
-    managers : dict
-        managers
-    nf: int
-        number of active flavor below threshold
-    q2: float
-        squared matching scale
-    is_backward: bool
-        True for backward matching
-    L: float
-        :math:`\ln(\mu_F^2 / m_h^2)`
-    is_msbar: bool
-        add the |MSbar| contribution
-    """
-
-    log_label = "Matching"
-    # complete list of possible matching operators labels
-    full_labels = [
-        *br.singlet_labels,
-        (br.matching_hplus_pid, 21),
-        (br.matching_hplus_pid, 100),
-        (21, br.matching_hplus_pid),
-        (100, br.matching_hplus_pid),
-        (br.matching_hplus_pid, br.matching_hplus_pid),
-        (200, 200),
-        (200, br.matching_hminus_pid),
-        (br.matching_hminus_pid, 200),
-        (br.matching_hminus_pid, br.matching_hminus_pid),
-    ]
-    # still valid in QED since Sdelta and Vdelta matchings are diagonal
-    full_labels_qed = copy.deepcopy(full_labels)
-
-    def __init__(self, config, managers, nf, q2, is_backward, L, is_msbar):
-        super().__init__(config, managers, Segment(q2, q2, nf))
-        self.backward_method = config["backward_inversion"] if is_backward else None
-        if is_backward:
-            self.is_intrinsic = True
-        else:
-            self.is_intrinsic = bool(len(config["intrinsic_range"]) != 0)
-        self.L = L
-        self.is_msbar = is_msbar
-        # Note for the moment only QCD matching is implemented
-        self.order = tuple(config["matching_order"])
-
-    @property
-    def labels(self):
-        """Necessary sector labels to compute.
-
-        Returns
-        -------
-        list(str)
-            sector labels
-        """
-        labels = []
-        # non-singlet labels
-        if self.config["debug_skip_non_singlet"]:
-            logger.warning("%s: skipping non-singlet sector", self.log_label)
-        else:
-            labels.append((200, 200))
-            if self.is_intrinsic or self.backward_method is not None:
-                # intrinsic labels, which are not zero at NLO
-                labels.append((br.matching_hminus_pid, br.matching_hminus_pid))
-                # These contributions are always 0 for the moment
-                # labels.extend([(200, br.matching_hminus_pid), (br.matching_hminus_pid, 200)])
-        # same for singlet
-        if self.config["debug_skip_singlet"]:
-            logger.warning("%s: skipping singlet sector", self.log_label)
-        else:
-            labels.extend(
-                [
-                    *br.singlet_labels,
-                    (br.matching_hplus_pid, 21),
-                    (br.matching_hplus_pid, 100),
-                ]
-            )
-            if self.is_intrinsic or self.backward_method is not None:
-                labels.extend(
-                    [
-                        (21, br.matching_hplus_pid),
-                        (100, br.matching_hplus_pid),
-                        (br.matching_hplus_pid, br.matching_hplus_pid),
-                    ]
-                )
-        return labels
-
-    # def quad_ker(self, label, logx, areas):
-    #     """Return partially initialized integrand function.
-
-    #     Parameters
-    #     ----------
-    #     label: tuple
-    #         operator element pids
-    #     logx: float
-    #         Mellin inversion point
-    #     areas : tuple
-    #         basis function configuration
-
-    #     Returns
-    #     -------
-    #     functools.partial
-    #         partially initialized integration kernel
-    #     """
-    #     return functools.partial(
-    #         quad_ker,
-    #         order=self.order,
-    #         mode0=label[0],
-    #         mode1=label[1],
-    #         is_log=self.int_disp.log,
-    #         logx=logx,
-    #         areas=areas,
-    #         a_s=self.a_s,
-    #         nf=self.nf,
-    #         L=self.L,
-    #         sv_mode=self.sv_mode,
-    #         Lsv=np.log(self.xif2),
-    #         backward_method=self.backward_method,
-    #         is_msbar=self.is_msbar,
-    #         is_polarized=self.config["polarized"],
-    #         is_time_like=self.config["time_like"],
-    #     )
-
-    def run_op_integration(self, log_grid):
-        """Run the integration for each grid point.
-
-        Parameters
-        ----------
-        log_grid : tuple(k, logx)
-            log grid point with relative index
-
-        Returns
-        -------
-        list
-            computed operators at the give grid point
-        """
-        column = []
-        k, logx = log_grid
-        # call(!) self.labels only once
-        labels = self.labels
-        start_time = time.perf_counter()
-        # start preparing C arguments
-        cfg = ekors.lib.empty_qcd_args()
-        cfg.order_qcd = self.order[0]
-        cfg.is_polarized = self.config["polarized"]
-        cfg.is_time_like = self.config["time_like"]
-        cfg.nf = self.nf
-        cfg.py = ekors.ffi.cast("void *", cb_quad_ker_ome.address)
-        cfg.is_log = self.int_disp.log
-        cfg.logx = logx
-        cfg.L = np.log(self.xif2)
-        # cfg.method_num = 1
-        cfg.as1 = self.as_list[1]
-        cfg.as0 = self.as_list[0]
-        # cfg.ev_op_iterations = self.config["ev_op_iterations"]
-        # cfg.ev_op_max_order_qcd = self.config["ev_op_max_order"][0]
-        # cfg.sv_mode_num = 1
-        # cfg.is_threshold = self.is_threshold
-
-        # iterate basis functions
-        for l, bf in enumerate(self.int_disp):
-            if k == l and l == self.grid_size - 1:
-                continue
-            # add emtpy labels with 0s
-            if bf.is_below_x(np.exp(logx)):
-                column.append({label: (0.0, 0.0) for label in labels})
-                continue
-            temp_dict = {}
-            # prepare areas for C
-            curareas = bf.areas_representation
-            areas_len = curareas.shape[0] * curareas.shape[1]
-            # force the variable in scope
-            areas_ffi = ekors.ffi.new(
-                f"double[{areas_len}]", curareas.flatten().tolist()
-            )
-            cfg.areas = areas_ffi
-            cfg.areas_x = curareas.shape[0]
-            cfg.areas_y = curareas.shape[1]
-            # iterate sectors
-            for label in labels:
-                cfg.mode0 = label[0]
-                cfg.mode1 = label[1]
-                # construct the low level object
-                func = LowLevelCallable(
-                    ekors.lib.rust_quad_ker_ome, ekors.ffi.addressof(cfg)
-                )
-                res = integrate.quad(
-                    func,
-                    0.5,
-                    1.0 - self._mellin_cut,
-                    epsabs=1e-12,
-                    epsrel=1e-5,
-                    limit=100,
-                    full_output=1,
-                )
-                temp_dict[label] = res[:2]
-            column.append(temp_dict)
-        logger.info(
-            "%s: computing operators - %u/%u took: %6f s",
-            self.log_label,
-            k + 1,
-            self.grid_size,
-            (time.perf_counter() - start_time),
-        )
-        return column
-
-    @property
-    def a_s(self):
-        """Return the computed values for :math:`a_s`.
-
-        Note that here you need to use :math:`a_s^{n_f+1}`
-        """
-        sc = self.managers["couplings"]
-        return sc.a_s(
-            self.q2_from
-            * (self.xif2 if self.sv_mode == sv.Modes.exponentiated else 1.0),
-            nf_to=self.nf + 1,
-        )
-
-    def compute(self):
-        """Compute the actual operators (i.e. run the integrations)."""
-        self.initialize_op_members()
-
-        # At LO you don't need anything else
-        if self.order[0] == 0:
-            logger.info("%s: no need to compute matching at LO", self.log_label)
-            return
-        logger.info(
-            "%s: order: (%d, %d), backward method: %s",
-            self.log_label,
-            self.order[0],
-            self.order[1],
-            self.backward_method,
-        )
-
-        self.integrate()
diff --git a/src/eko/evolution_operator/quad_ker.py b/src/eko/evolution_operator/quad_ker.py
index 798cebe6b..94772ef72 100644
--- a/src/eko/evolution_operator/quad_ker.py
+++ b/src/eko/evolution_operator/quad_ker.py
@@ -13,6 +13,7 @@
 from ..kernels import singlet as s
 from ..matchings import Segment
 from ..member import OpMember
+from .operator_matrix_element import build_ome
 
 logger = logging.getLogger(__name__)
 
@@ -66,6 +67,7 @@ def select_singlet_element(ker, mode0, mode1):
         nb.types.uint,  # ev_op_max_order_qcd
         nb.types.uint,  # sv_mode_num
         nb.types.bool_,  # is_threshold
+        nb.types.double,  # Lsv
     ),
     cache=True,
     nopython=True,
@@ -95,6 +97,7 @@ def cb_quad_ker_qcd(
     ev_op_max_order_qcd,
     _sv_mode_num,
     is_threshold,
+    Lsv,  # dummy variable
 ):
     """C Callback inside integration kernel."""
     # recover complex variables
@@ -158,6 +161,105 @@ def cb_quad_ker_qcd(
     return np.real(res)
 
 
+@nb.cfunc(
+    nb.types.double(
+        nb.types.CPointer(nb.types.double),  # re_gamma_ns_raw
+        nb.types.CPointer(nb.types.double),  # im_gamma_ns_raw
+        nb.types.double,  # re_n
+        nb.types.double,  # im_n
+        nb.types.double,  # re_jac
+        nb.types.double,  # im_jac
+        nb.types.uint,  # order_qcd
+        nb.types.bool_,  # is_singlet
+        nb.types.uint,  # mode0
+        nb.types.uint,  # mode1
+        nb.types.uint,  # nf
+        nb.types.bool_,  # is_log
+        nb.types.double,  # logx
+        nb.types.CPointer(nb.types.double),  # areas_raw
+        nb.types.uint,  # areas_x
+        nb.types.uint,  # areas_y
+        nb.types.double,  # L
+        nb.types.uint,  # method_num
+        nb.types.double,  # as1
+        nb.types.double,  # as0
+        nb.types.uint,  # ev_op_iterations
+        nb.types.uint,  # ev_op_max_order_qcd
+        nb.types.uint,  # sv_mode_num
+        nb.types.bool_,  # is_threshold
+        nb.types.double,  # Lsv
+    ),
+    cache=True,
+    nopython=True,
+)
+def cb_quad_ker_ome(
+    re_gamma_raw,
+    im_gamma_raw,
+    re_n,
+    im_n,
+    re_jac,
+    im_jac,
+    order_qcd,
+    is_singlet,
+    mode0,
+    mode1,
+    nf,
+    is_log,
+    logx,
+    areas_raw,
+    areas_x,
+    areas_y,
+    L,
+    _method_num,  # dummy variable
+    as1,
+    as0,
+    ev_op_iterations,  # dummy variable
+    ev_op_max_order_qcd,  # dummy variable
+    _sv_mode_num,  # dummy variable
+    is_threshold,  # dummy variable
+    Lsv,
+):
+    """C Callback inside integration kernel."""
+    # recover complex variables
+    n = re_n + im_n * 1j
+    jac = re_jac + im_jac * 1j
+    # compute basis functions
+    areas = nb.carray(areas_raw, (areas_x, areas_y))
+    pj = interpolation.evaluate_grid(n, is_log, logx, areas)
+    # TODO recover parameters
+    sv_mode = sv.Modes.exponentiated
+    order = (order_qcd, 0)
+    if is_singlet:
+        indices = {21: 0, 100: 1, 90: 2}
+        # reconstruct singlet matrices
+        re_ome_singlet = nb.carray(re_ome_raw, (order_qcd, 3, 3))
+        im_ome_singlet = nb.carray(im_ome_raw, (order_qcd, 3, 3))
+        A = re_ome_singlet + im_ome_singlet * 1j
+    else:
+        indices = {200: 0, 91: 1}
+        # construct non-singlet matrices
+        re_ome_ns = nb.carray(re_ome_raw, (order_qcd, 2, 2))
+        im_ome_ns = nb.carray(im_ome_raw, (order_qcd, 2, 2, 2))
+        A = re_ome_ns + im_ome_ns * 1j
+
+    # correct for scale variations
+    if sv_mode == sv.Modes.exponentiated:
+        A = sv.exponentiated.gamma_variation(A, order, nf, Lsv)
+
+    # TODO recover InversionMethod
+    backward_method = "exact"
+
+    # build the expansion in alpha_s depending on the strategy
+    ker = build_ome(A, order, a_s, backward_method)
+
+    # select the needed matrix element
+    ker = ker[indices[mode0], indices[mode1]]
+
+    # recombine everything
+    res = ker * pj * jac
+    return np.real(res)
+
+
 # from ..kernels import singlet_qed as qed_s
 # from ..kernels import non_singlet_qed as qed_ns
 # from ..kernels import valence_qed as qed_v
diff --git a/src/eko/evolution_operator/quad_ker_ome.py b/src/eko/evolution_operator/quad_ker_ome.py
deleted file mode 100644
index 5b0a9378a..000000000
--- a/src/eko/evolution_operator/quad_ker_ome.py
+++ /dev/null
@@ -1,169 +0,0 @@
-r"""Integration kernels."""
-
-import logging
-
-import numba as nb
-import numpy as np
-from scipy import LowLevelCallable, integrate
-
-from .. import basis_rotation as br
-from .. import interpolation
-from .. import scale_variations as sv
-from ..io.types import InversionMethod
-from ..kernels import non_singlet as ns
-from ..kernels import singlet as s
-from ..matchings import Segment
-from ..member import OpMember
-
-logger = logging.getLogger(__name__)
-
-
-@nb.njit(cache=True)
-def build_ome(A, matching_order, a_s, backward_method):
-    r"""Construct the matching expansion in :math:`a_s` with the appropriate method.
-
-    Parameters
-    ----------
-    A : numpy.ndarray
-        list of |OME|
-    matching_order : tuple(int,int)
-        perturbation matching order
-    a_s : float
-        strong coupling, needed only for the exact inverse
-    backward_method : InversionMethod or None
-        empty or method for inverting the matching condition (exact or expanded)
-
-    Returns
-    -------
-    ome : numpy.ndarray
-        matching operator matrix
-
-    """
-    # to get the inverse one can use this FORM snippet
-    # Symbol a;
-    # NTensor c,d,e;
-    # Local x=-(a*c+a**2* d + a**3 * e);
-    # Local bi = 1+x+x**2+x**3;
-    # Print;
-    # .end
-    ome = np.eye(len(A[0]), dtype=np.complex_)
-    A = A[:, :, :]
-    A = np.ascontiguousarray(A)
-
-    if backward_method is InversionMethod.EXPANDED:
-        # expended inverse
-        if matching_order[0] >= 1:
-            ome -= a_s * A[0]
-        if matching_order[0] >= 2:
-            ome += a_s**2 * (-A[1] + A[0] @ A[0])
-        if matching_order[0] >= 3:
-            ome += a_s**3 * (-A[2] + A[0] @ A[1] + A[1] @ A[0] - A[0] @ A[0] @ A[0])
-    else:
-        # forward or exact inverse
-        if matching_order[0] >= 1:
-            ome += a_s * A[0]
-        if matching_order[0] >= 2:
-            ome += a_s**2 * A[1]
-        if matching_order[0] >= 3:
-            ome += a_s**3 * A[2]
-        # need inverse exact ?  so add the missing pieces
-        if backward_method is InversionMethod.EXACT:
-            ome = np.linalg.inv(ome)
-    return ome
-
-
-@nb.cfunc(
-    nb.types.double(
-        nb.types.CPointer(nb.types.double),  # re_gamma_ns_raw
-        nb.types.CPointer(nb.types.double),  # im_gamma_ns_raw
-        nb.types.double,  # re_n
-        nb.types.double,  # im_n
-        nb.types.double,  # re_jac
-        nb.types.double,  # im_jac
-        nb.types.uint,  # order_qcd
-        nb.types.bool_,  # is_singlet
-        nb.types.uint,  # mode0
-        nb.types.uint,  # mode1
-        nb.types.uint,  # nf
-        nb.types.bool_,  # is_log
-        nb.types.double,  # logx
-        nb.types.CPointer(nb.types.double),  # areas_raw
-        nb.types.uint,  # areas_x
-        nb.types.uint,  # areas_y
-        nb.types.double,  # L
-        nb.types.uint,  # method_num
-        nb.types.double,  # as1
-        nb.types.double,  # as0
-        nb.types.uint,  # ev_op_iterations
-        nb.types.uint,  # ev_op_max_order_qcd
-        nb.types.uint,  # sv_mode_num
-        nb.types.bool_,  # is_threshold
-    ),
-    cache=True,
-    nopython=True,
-)
-def cb_quad_ker_ome(
-    re_gamma_raw,
-    im_gamma_raw,
-    re_n,
-    im_n,
-    re_jac,
-    im_jac,
-    order_qcd,
-    is_singlet,
-    mode0,
-    mode1,
-    nf,
-    is_log,
-    logx,
-    areas_raw,
-    areas_x,
-    areas_y,
-    L,
-    _method_num,  # dummy variable
-    as1,
-    as0,
-    ev_op_iterations,  # dummy variable
-    ev_op_max_order_qcd,  # dummy variable
-    _sv_mode_num,  # dummy variable
-    is_threshold,  # dummy variable
-):
-    """C Callback inside integration kernel."""
-    # recover complex variables
-    n = re_n + im_n * 1j
-    jac = re_jac + im_jac * 1j
-    # compute basis functions
-    areas = nb.carray(areas_raw, (areas_x, areas_y))
-    pj = interpolation.evaluate_grid(n, is_log, logx, areas)
-    # TODO recover parameters
-    sv_mode = sv.Modes.exponentiated
-    order = (order_qcd, 0)
-    if is_singlet:
-        indices = {21: 0, 100: 1, 90: 2}
-        # reconstruct singlet matrices
-        re_ome_singlet = nb.carray(re_ome_raw, (order_qcd, 3, 3))
-        im_ome_singlet = nb.carray(im_ome_raw, (order_qcd, 3, 3))
-        A = re_ome_singlet + im_ome_singlet * 1j
-    else:
-        indices = {200: 0, 91: 1}
-        # construct non-singlet matrices
-        re_ome_ns = nb.carray(re_ome_raw, (order_qcd, 2, 2))
-        im_ome_ns = nb.carray(im_ome_raw, (order_qcd, 2, 2, 2))
-        A = re_ome_ns + im_ome_ns * 1j
-
-    # correct for scale variations
-    if sv_mode == sv.Modes.exponentiated:
-        A = sv.exponentiated.gamma_variation(A, order, nf, L)
-
-    # TODO recover InversionMethod
-    backward_method = "exact"
-
-    # build the expansion in alpha_s depending on the strategy
-    ker = build_ome(A, order, a_s, backward_method)
-
-    # select the needed matrix element
-    ker = ker[indices[mode0], indices[mode1]]
-
-    # recombine everything
-    res = ker * pj * jac
-    return np.real(res)