diff --git a/benchmarks/eko/benchmark_evol_to_unity.py b/benchmarks/eko/benchmark_evol_to_unity.py
index 6baf6bfa8..978d11c62 100644
--- a/benchmarks/eko/benchmark_evol_to_unity.py
+++ b/benchmarks/eko/benchmark_evol_to_unity.py
@@ -11,7 +11,7 @@
 from eko.io.runcards import OperatorCard, TheoryCard
 from eko.runner.legacy import Runner
 
-# from eko.matching_conditions.operator_matrix_element import OperatorMatrixElement
+# from ekore.matching_conditions.operator_matrix_element import OperatorMatrixElement
 
 
 def update_cards(theory: TheoryCard, operator: OperatorCard):
diff --git a/benchmarks/eko/benchmark_ad.py b/benchmarks/ekore/benchmark_ad.py
similarity index 98%
rename from benchmarks/eko/benchmark_ad.py
rename to benchmarks/ekore/benchmark_ad.py
index cb5c07e66..08919150e 100644
--- a/benchmarks/eko/benchmark_ad.py
+++ b/benchmarks/ekore/benchmark_ad.py
@@ -2,8 +2,8 @@
 import numpy as np
 import pytest
 
-import eko.anomalous_dimensions.as2 as ad_as2
-import eko.harmonics as h
+import ekore.anomalous_dimensions.unpolarized.space_like.as2 as ad_as2
+import ekore.harmonics as h
 from eko.constants import CA, CF, TR
 
 
diff --git a/benchmarks/performance/harmonics.py b/benchmarks/performance/harmonics.py
index f803bcda0..c0b7c4a95 100644
--- a/benchmarks/performance/harmonics.py
+++ b/benchmarks/performance/harmonics.py
@@ -1,6 +1,6 @@
 import numpy as np
 
-from eko import anomalous_dimensions as ad
+import ekore.anomalous_dimensions.unpolarized.space_like as ad
 from eko.mellin import Path
 
 NF = 5
diff --git a/doc/source/overview/features.rst b/doc/source/overview/features.rst
index 719bc5d61..3edf80aff 100644
--- a/doc/source/overview/features.rst
+++ b/doc/source/overview/features.rst
@@ -1,8 +1,8 @@
 Features
 ========
 
-- perturbation orders: :mod:`LO <eko.anomalous_dimensions.lo>` + :mod:`NLO <eko.anomalous_dimensions.nlo>`
-  + :mod:`NNLO <eko.anomalous_dimensions.nnlo>`
+- perturbation orders: :mod:`LO <ekore.anomalous_dimensions.lo>` + :mod:`NLO <ekore.anomalous_dimensions.nlo>`
+  + :mod:`NNLO <ekore.anomalous_dimensions.nnlo>`
 - evolution in an (almost) arbitrary sequence of |FNS| (:class:`~eko.thresholds.ThresholdsAtlas`)
 - :ref:`theory/pQCD:scale variations`
 - different :doc:`solutions </theory/DGLAP>` of |DGLAP| differential equation
diff --git a/doc/source/theory/Mellin.rst b/doc/source/theory/Mellin.rst
index 083f3e9a9..2c160902f 100644
--- a/doc/source/theory/Mellin.rst
+++ b/doc/source/theory/Mellin.rst
@@ -124,7 +124,7 @@ the Mellin inverse.
     S_m(N) = \sum\limits_{j=1}^N \frac 1 {j^m} = \frac{(-1)^{m-1}}{(m-1)!} \psi_{m-1}(N+1)+c_m \quad
     \text{with},\quad c_m = \left\{\begin{array}{ll} \gamma_E, & m=1\\ \zeta(m), & m>1\end{array} \right.
 
-  where :math:`\psi_k(N)` is the :math:`k`-th polygamma function (implemented as :meth:`~eko.harmonics.polygamma.cern_polygamma`)
+  where :math:`\psi_k(N)` is the :math:`k`-th polygamma function (implemented as :meth:`~ekore.harmonics.polygamma.cern_polygamma`)
   and :math:`\zeta` the Riemann zeta function (using :func:`scipy.special.zeta`).
 
 - for the sums :math:`S_{-m}(N)` and m > 0 we use :cite:`Gluck:1989ze`:
@@ -149,7 +149,7 @@ the Mellin inverse.
   which express higher weight sums in terms of simple one :math:`S_{m}, S_{-m}` and some irreducible integrals.
   The above prescription on the analytical continuation of :math:`\eta` is applied.
 
-The complete list of harmonics sums available in :mod:`eko.harmonics` is:
+The complete list of harmonics sums available in :mod:`ekore.harmonics` is:
 
     - weight 1:
 
@@ -166,20 +166,20 @@ The complete list of harmonics sums available in :mod:`eko.harmonics` is:
         .. math::
             S_{3}, S_{2,1}, S_{2,-1}, S_{-2,1}, S_{-2,-1}, S_{-3}
 
-        these sums relies on the integrals :mod:`eko.harmonics.g_functions` :cite:`MuselliPhD,Blumlein:1998if`
+        these sums relies on the integrals :mod:`ekore.harmonics.g_functions` :cite:`MuselliPhD,Blumlein:1998if`
 
     - weight 4:
 
         .. math ::
             S_{4}, S_{3,1}, S_{2,1,1}, S_{-2,-2}, S_{-3, 1}, S_{-4}
 
-        these sums relies on the integrals :mod:`eko.harmonics.g_functions` :cite:`MuselliPhD,Blumlein:1998if`
+        these sums relies on the integrals :mod:`ekore.harmonics.g_functions` :cite:`MuselliPhD,Blumlein:1998if`
 
     - weight 5:
 
         .. math ::
             S_{5}, S_{4,1}, S_{3,1,1}, S_{2,3}, S_{2,2,1}, S_{2,1,1,1}, S_{2,1,-2}, S_{2,-3}, S_{-2,3}, S_{-2,2,1}, S_{-2,1,1,1}, S_{-5}
 
-        these sums relies on the integrals :mod:`eko.harmonics.f_functions` :cite:`Blumlein:2009ta`
+        these sums relies on the integrals :mod:`ekore.harmonics.f_functions` :cite:`Blumlein:2009ta`
 
-We have also implemented a recursive computation of simple harmonics (single index), see :func:`eko.harmonics.polygamma.recursive_harmonic_sum`
+We have also implemented a recursive computation of simple harmonics (single index), see :func:`ekore.harmonics.polygamma.recursive_harmonic_sum`
diff --git a/pyproject.toml b/pyproject.toml
index f290ce339..79889abbe 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -133,6 +133,7 @@ python_functions = ['test_*', 'benchmark_*']
 addopts = [
   '--cov=eko',
   '--cov=ekobox',
+  '--cov=ekore',
   '--cov-report=html',
   '--cov-report=xml',
   '--strict-markers',
diff --git a/src/eko/beta.py b/src/eko/beta.py
index f83d013bb..92212bf58 100644
--- a/src/eko/beta.py
+++ b/src/eko/beta.py
@@ -7,7 +7,7 @@
 import numba as nb
 
 from . import constants
-from .harmonics.constants import zeta3
+from ekore.harmonics.constants import zeta3
 
 
 @nb.njit(cache=True)
diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index 5818c90fc..ee17b9cd6 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -13,7 +13,10 @@
 import numpy as np
 from scipy import integrate
 
-from .. import anomalous_dimensions as ad
+import ekore.anomalous_dimensions.polarized.space_like as ad_ps
+import ekore.anomalous_dimensions.unpolarized.space_like as ad_us
+import ekore.anomalous_dimensions.unpolarized.time_like as ad_ut
+
 from .. import basis_rotation as br
 from .. import interpolation, mellin
 from .. import scale_variations as sv
@@ -129,6 +132,8 @@ def quad_ker(
     ev_op_max_order,
     sv_mode,
     is_threshold,
+    is_polarized,
+    is_time_like,
 ):
     """Raw evolution kernel inside quad.
 
@@ -165,7 +170,11 @@ def quad_ker(
     sv_mode: int, `enum.IntEnum`
         scale variation mode, see `eko.scale_variations.Modes`
     is_threshold : boolean
-        is this an itermediate threshold operator?
+        is this an intermediate threshold operator?
+    is_polarized : boolean
+        is polarized evolution ?
+    is_time_like : boolean
+        is time-like evolution ?
 
     Returns
     -------
@@ -179,7 +188,16 @@ def quad_ker(
 
     # compute the actual evolution kernel
     if ker_base.is_singlet:
-        gamma_singlet = ad.gamma_singlet(order, ker_base.n, nf)
+        if is_polarized:
+            if is_time_like:
+                raise NotImplementedError("Polarized, time-like is not implemented")
+            else:
+                gamma_singlet = ad_ps.gamma_singlet(order, ker_base.n, nf)
+        else:
+            if is_time_like:
+                gamma_singlet = ad_ut.gamma_singlet(order, ker_base.n, nf)
+            else:
+                gamma_singlet = ad_us.gamma_singlet(order, ker_base.n, nf)
         # scale var exponentiated is directly applied on gamma
         if sv_mode == sv.Modes.exponentiated:
             gamma_singlet = sv.exponentiated.gamma_variation(
@@ -202,7 +220,16 @@ def quad_ker(
             ) @ np.ascontiguousarray(ker)
         ker = select_singlet_element(ker, mode0, mode1)
     else:
-        gamma_ns = ad.gamma_ns(order, mode0, ker_base.n, nf)
+        if is_polarized:
+            if is_time_like:
+                raise NotImplementedError("Polarized, time-like is not implemented")
+            else:
+                gamma_ns = ad_ps.gamma_ns(order, mode0, ker_base.n, nf)
+        else:
+            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)
         if sv_mode == sv.Modes.exponentiated:
             gamma_ns = sv.exponentiated.gamma_variation(gamma_ns, order, nf, L)
         ker = ns.dispatcher(
@@ -381,6 +408,8 @@ def quad_ker(self, label, logx, areas):
             ev_op_max_order=tuple(self.config["ev_op_max_order"]),
             sv_mode=self.sv_mode,
             is_threshold=self.is_threshold,
+            is_polarized=self.config["polarized"],
+            is_time_like=self.config["time_like"],
         )
 
     def initialize_op_members(self):
diff --git a/src/eko/evolution_operator/grid.py b/src/eko/evolution_operator/grid.py
index 1b55a2072..270d434a8 100644
--- a/src/eko/evolution_operator/grid.py
+++ b/src/eko/evolution_operator/grid.py
@@ -11,12 +11,12 @@
 import numpy as np
 import numpy.typing as npt
 
-from .. import matching_conditions, member
+from .. import member
 from .. import scale_variations as sv
 from ..io.runcards import Configs, Debug
-from ..matching_conditions.operator_matrix_element import OperatorMatrixElement
 from ..thresholds import flavor_shift, is_downward_path
-from . import Operator, flavors, physical
+from . import Operator, flavors, matching_condition, physical
+from .operator_matrix_element import OperatorMatrixElement
 
 logger = logging.getLogger(__name__)
 
@@ -89,6 +89,8 @@ def __init__(
         config["n_integration_cores"] = configs.n_integration_cores
         config["debug_skip_singlet"] = debug.skip_singlet
         config["debug_skip_non_singlet"] = debug.skip_non_singlet
+        config["polarized"] = configs.polarized
+        config["time_like"] = configs.time_like
 
         if method not in [
             "iterate-exact",
@@ -244,7 +246,7 @@ def generate(self, q2):
 
             # join with the basis rotation, since matching requires c+ (or likewise)
             if is_downward:
-                matching = matching_conditions.MatchingCondition.split_ad_to_evol_map(
+                matching = matching_condition.MatchingCondition.split_ad_to_evol_map(
                     self._matching_operators[op.q2_to],
                     op.nf - 1,
                     op.q2_to,
@@ -255,7 +257,7 @@ def generate(self, q2):
                 )
                 final_op = final_op @ matching @ invrot @ phys_op
             else:
-                matching = matching_conditions.MatchingCondition.split_ad_to_evol_map(
+                matching = matching_condition.MatchingCondition.split_ad_to_evol_map(
                     self._matching_operators[op.q2_to],
                     op.nf,
                     op.q2_to,
diff --git a/src/eko/matching_conditions/__init__.py b/src/eko/evolution_operator/matching_condition.py
similarity index 100%
rename from src/eko/matching_conditions/__init__.py
rename to src/eko/evolution_operator/matching_condition.py
diff --git a/src/eko/matching_conditions/operator_matrix_element.py b/src/eko/evolution_operator/operator_matrix_element.py
similarity index 62%
rename from src/eko/matching_conditions/operator_matrix_element.py
rename to src/eko/evolution_operator/operator_matrix_element.py
index ed0299294..a1c7c445b 100644
--- a/src/eko/matching_conditions/operator_matrix_element.py
+++ b/src/eko/evolution_operator/operator_matrix_element.py
@@ -1,7 +1,4 @@
-"""
-This module defines the |OME| for the non-trivial matching conditions in the
-|VFNS| evolution.
-"""
+"""The |OME| for the non-trivial matching conditions in the |VFNS| evolution."""
 
 import functools
 import logging
@@ -9,116 +6,36 @@
 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 ekore import harmonics
+
 from .. import basis_rotation as br
-from .. import harmonics
-from ..evolution_operator import Operator, QuadKerBase
-from . import as1, as2, as3
+from . import Operator, QuadKerBase
 
 logger = logging.getLogger(__name__)
 
 
-@nb.njit(cache=True)
-def A_singlet(matching_order, n, sx, nf, L, is_msbar, sx_ns=None):
-    r"""
-    Computes the tower of the singlet |OME|.
-
-    Parameters
-    ----------
-        matching_order : tuple(int,int)
-            perturbative matching_order
-        n : complex
-            Mellin variable
-        sx : list
-            singlet like harmonic sums cache
-        nf: int
-            number of active flavor below threshold
-        L : float
-            :math:``\ln(\mu_F^2 / m_h^2)``
-        is_msbar: bool
-            add the |MSbar| contribution
-        sx_ns : list
-            non-singlet like harmonic sums cache
-
-    Returns
-    -------
-        A_singlet : numpy.ndarray
-            singlet |OME|
-
-    See Also
-    --------
-        eko.matching_conditions.nlo.A_singlet_1 : :math:`A^{S,(1)}(N)`
-        eko.matching_conditions.nlo.A_hh_1 : :math:`A_{HH}^{(1)}(N)`
-        eko.matching_conditions.nlo.A_gh_1 : :math:`A_{gH}^{(1)}(N)`
-        eko.matching_conditions.nnlo.A_singlet_2 : :math:`A_{S,(2)}(N)`
-    """
-    A_s = np.zeros((matching_order[0], 3, 3), np.complex_)
-    if matching_order[0] >= 1:
-        A_s[0] = as1.A_singlet(n, sx, L)
-    if matching_order[0] >= 2:
-        A_s[1] = as2.A_singlet(n, sx, L, is_msbar)
-    if matching_order[0] >= 3:
-        A_s[2] = as3.A_singlet(n, sx, sx_ns, nf, L)
-    return A_s
-
-
-@nb.njit(cache=True)
-def A_non_singlet(matching_order, n, sx, nf, L):
-    r"""
-    Computes the tower of the non-singlet |OME|
-
-    Parameters
-    ----------
-        matching_order : tuple(int,int)
-            perturbative matching_order
-        n : complex
-            Mellin variable
-        sx : list
-            harmonic sums cache
-        nf: int
-            number of active flavor below threshold
-        L : float
-            :math:``\ln(\mu_F^2 / m_h^2)``
-
-    Returns
-    -------
-        A_non_singlet : numpy.ndarray
-            non-singlet |OME|
-
-    See Also
-    --------
-        eko.matching_conditions.nlo.A_hh_1 : :math:`A_{HH}^{(1)}(N)`
-        eko.matching_conditions.nnlo.A_ns_2 : :math:`A_{qq,H}^{NS,(2)}`
-    """
-    A_ns = np.zeros((matching_order[0], 2, 2), np.complex_)
-    if matching_order[0] >= 1:
-        A_ns[0] = as1.A_ns(n, sx, L)
-    if matching_order[0] >= 2:
-        A_ns[1] = as2.A_ns(n, sx, L)
-    if matching_order[0] >= 3:
-        A_ns[2] = as3.A_ns(n, sx, nf, L)
-    return A_ns
-
-
 @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.
+    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 : ["exact", "expanded" or ""]
-            empty or method for inverting the matching condition (exact or expanded)
+    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 : ["exact", "expanded" or ""]
+        empty or method for inverting the matching condition (exact or expanded)
 
     Returns
     -------
-        ome : numpy.ndarray
-            matching operator matrix
+    ome : numpy.ndarray
+        matching operator matrix
     """
     # to get the inverse one can use this FORM snippet
     # Symbol a;
@@ -153,41 +70,58 @@ def build_ome(A, matching_order, a_s, backward_method):
 
 @nb.njit(cache=True)
 def quad_ker(
-    u, order, mode0, mode1, is_log, logx, areas, a_s, nf, L, backward_method, is_msbar
+    u,
+    order,
+    mode0,
+    mode1,
+    is_log,
+    logx,
+    areas,
+    a_s,
+    nf,
+    L,
+    backward_method,
+    is_msbar,
+    is_polarized,
+    is_time_like,
 ):
-    r"""
-    Raw kernel inside quad
+    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 : ["exact", "expanded" or ""]
-            empty or method for inverting the matching condition (exact or expanded)
-        is_msbar: bool
-            add the |MSbar| contribution
+    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 : ["exact", "expanded" or ""]
+        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 : float
+        evaluated integration kernel
     """
     ker_base = QuadKerBase(u, is_log, logx, mode0)
     integrand = ker_base.integrand(areas)
@@ -220,10 +154,28 @@ def quad_ker(
     # compute the ome
     if ker_base.is_singlet:
         indices = {21: 0, 100: 1, 90: 2}
-        A = A_singlet(order, ker_base.n, sx, nf, L, is_msbar, sx_ns)
+        if is_polarized:
+            if is_time_like:
+                raise NotImplementedError("Polarized, time-like is not implemented")
+            else:
+                A = ome_ps.A_singlet(order, ker_base.n, sx, nf, L, is_msbar, sx_ns)
+        else:
+            if is_time_like:
+                A = ome_ut.A_singlet(order, ker_base.n, sx, nf, L, is_msbar, sx_ns)
+            else:
+                A = ome_us.A_singlet(order, ker_base.n, sx, nf, L, is_msbar, sx_ns)
     else:
         indices = {200: 0, 91: 1}
-        A = A_non_singlet(order, ker_base.n, sx, nf, L)
+        if is_polarized:
+            if is_time_like:
+                raise NotImplementedError("Polarized, time-like is not implemented")
+            else:
+                A = ome_us.A_non_singlet(order, ker_base.n, sx, nf, L)
+        else:
+            if is_time_like:
+                A = ome_us.A_non_singlet(order, ker_base.n, sx, nf, L)
+            else:
+                A = ome_us.A_non_singlet(order, ker_base.n, sx, nf, L)
 
     # build the expansion in alpha_s depending on the strategy
     ker = build_ome(A, order, a_s, backward_method)
@@ -287,14 +239,13 @@ def __init__(self, config, managers, nf, q2, is_backward, L, is_msbar):
 
     @property
     def labels(self):
-        """Computes the necessary sector labels to compute.
+        """Necessary sector labels to compute.
 
         Returns
         -------
         list(str)
             sector labels
         """
-
         labels = []
         # non-singlet labels
         if self.config["debug_skip_non_singlet"]:
@@ -328,7 +279,7 @@ def labels(self):
         return labels
 
     def quad_ker(self, label, logx, areas):
-        """Partially initialized integrand function.
+        """Return partially initialized integrand function.
 
         Parameters
         ----------
@@ -357,11 +308,13 @@ def quad_ker(self, label, logx, areas):
             L=self.L,
             backward_method=self.backward_method,
             is_msbar=self.is_msbar,
+            is_polarized=self.config["polarized"],
+            is_time_like=self.config["time_like"],
         )
 
     @property
     def a_s(self):
-        """Returns the computed values for :math:`a_s`.
+        """Compute values for :math:`a_s`.
 
         Note that here you need to use :math:`a_s^{n_f+1}`
         """
@@ -371,7 +324,7 @@ def a_s(self):
         )
 
     def compute(self):
-        """Compute the actual operators (i.e. run the integrations)"""
+        """Compute the actual operators (i.e. run the integrations)."""
         self.initialize_op_members()
 
         # At LO you don't need anything else
diff --git a/src/eko/gamma.py b/src/eko/gamma.py
index cd4e9aa3e..b2c4a857b 100644
--- a/src/eko/gamma.py
+++ b/src/eko/gamma.py
@@ -5,7 +5,7 @@
 """
 import numba as nb
 
-from .harmonics.constants import zeta3, zeta4, zeta5
+from ekore.harmonics.constants import zeta3, zeta4, zeta5
 
 
 @nb.njit(cache=True)
diff --git a/src/eko/io/runcards.py b/src/eko/io/runcards.py
index 23cc20fa2..2ab9a7faa 100644
--- a/src/eko/io/runcards.py
+++ b/src/eko/io/runcards.py
@@ -88,6 +88,10 @@ class Configs(DictLike):
     r"""Whether to use polynomials in :math:`\log(x)`.
     If `false`, polynomials are in :math:`x`.
     """
+    polarized: bool
+    """If `true` do polarized evolution."""
+    time_like: bool
+    """If `true` do time-like evolution."""
     scvar_method: Optional[ScaleVariationsMethod]
     """"""
     inversion_method: Optional[InversionMethod]
@@ -326,6 +330,8 @@ def new_operator(self):
             "interpolation_is_log",
             "ev_op_iterations",
             "n_integration_cores",
+            "polarized",
+            "time_like",
         ):
             new["configs"][k] = old[k]
         max_order = old["ev_op_max_order"]
diff --git a/src/eko/kernels/singlet.py b/src/eko/kernels/singlet.py
index ab3d4dfc8..2bab5d35c 100644
--- a/src/eko/kernels/singlet.py
+++ b/src/eko/kernels/singlet.py
@@ -3,7 +3,7 @@
 import numba as nb
 import numpy as np
 
-from .. import anomalous_dimensions as ad
+from ekore import anomalous_dimensions as ad
 from .. import beta
 from . import as4_evolution_integrals as as4_ei
 from . import evolution_integrals as ei
diff --git a/src/ekobox/cards.py b/src/ekobox/cards.py
index 52569909b..231dc7266 100644
--- a/src/ekobox/cards.py
+++ b/src/ekobox/cards.py
@@ -38,6 +38,8 @@
         scvar_method=None,
         inversion_method=None,
         n_integration_cores=0,
+        polarized=False,
+        time_like=False,
     ),
     debug=dict(
         skip_singlet=False,
diff --git a/src/ekomark/data/db.py b/src/ekomark/data/db.py
index 782dbef30..a9fe0fd9f 100644
--- a/src/ekomark/data/db.py
+++ b/src/ekomark/data/db.py
@@ -22,3 +22,5 @@ class Operator(Base):  # pylint: disable=too-few-public-methods
     ev_op_iterations = Column(Integer)
     Q2grid = Column(Text)
     backward_inversion = Column(Text)
+    polarized = Column(Boolean)
+    time_like = Column(Boolean)
diff --git a/src/ekomark/data/operators.py b/src/ekomark/data/operators.py
index 4a161c0fd..bf19cb465 100644
--- a/src/ekomark/data/operators.py
+++ b/src/ekomark/data/operators.py
@@ -22,6 +22,8 @@
             targetgrid=None,
             inputpids=None,
             targetpids=None,
+            polarized=False,
+            time_like=False,
         ).items()
     )
 )
diff --git a/src/ekomark/navigator/navigator.py b/src/ekomark/navigator/navigator.py
index 48148de16..137ecca64 100644
--- a/src/ekomark/navigator/navigator.py
+++ b/src/ekomark/navigator/navigator.py
@@ -81,6 +81,8 @@ def fill_operators(self, op, obj):
         obj["iters"] = op["ev_op_iterations"]
         obj["skip_ns"] = op["debug_skip_non_singlet"]
         obj["skip_s"] = op["debug_skip_singlet"]
+        obj["pol"] = op["polarized"]
+        obj["time"] = op["time_like"]
 
     def fill_cache(self, cac, obj):
         """
diff --git a/tests/eko/anomalous_dimensions/__init__.py b/src/ekore/__init__.py
similarity index 100%
rename from tests/eko/anomalous_dimensions/__init__.py
rename to src/ekore/__init__.py
diff --git a/src/ekore/anomalous_dimensions/__init__.py b/src/ekore/anomalous_dimensions/__init__.py
new file mode 100644
index 000000000..47a877c03
--- /dev/null
+++ b/src/ekore/anomalous_dimensions/__init__.py
@@ -0,0 +1,62 @@
+r"""
+This module contains the Altarelli-Parisi splitting kernels.
+
+Normalization is given by
+
+.. math::
+    \mathbf{P}(x) = \sum\limits_{j=0} a_s^{j+1} \mathbf P^{(j)}(x)
+
+with :math:`a_s = \frac{\alpha_S(\mu^2)}{4\pi}`.
+The 3-loop references for the non-singlet :cite:`Moch:2004pa`
+and singlet :cite:`Vogt:2004mw` case contain also the lower
+order results. The results are also determined in Mellin space in
+terms of the anomalous dimensions (note the additional sign!)
+
+.. math::
+    \gamma(N) = - \mathcal{M}[\mathbf{P}(x)](N)
+"""
+
+import numba as nb
+import numpy as np
+
+
+@nb.njit(cache=True)
+def exp_singlet(gamma_S):
+    r"""
+    Computes the exponential and the eigensystem of the singlet anomalous dimension matrix
+
+    Parameters
+    ----------
+        gamma_S : numpy.ndarray
+            singlet anomalous dimension matrix
+
+    Returns
+    -------
+        exp : numpy.ndarray
+            exponential of the singlet anomalous dimension matrix :math:`\gamma_{S}(N)`
+        lambda_p : complex
+            positive eigenvalue of the singlet anomalous dimension matrix
+            :math:`\gamma_{S}(N)`
+        lambda_m : complex
+            negative eigenvalue of the singlet anomalous dimension matrix
+            :math:`\gamma_{S}(N)`
+        e_p : numpy.ndarray
+            projector for the positive eigenvalue of the singlet anomalous
+            dimension matrix :math:`\gamma_{S}(N)`
+        e_m : numpy.ndarray
+            projector for the negative eigenvalue of the singlet anomalous
+            dimension matrix :math:`\gamma_{S}(N)`
+    """
+    # compute eigenvalues
+    det = np.sqrt(
+        np.power(gamma_S[0, 0] - gamma_S[1, 1], 2) + 4.0 * gamma_S[0, 1] * gamma_S[1, 0]
+    )
+    lambda_p = 1.0 / 2.0 * (gamma_S[0, 0] + gamma_S[1, 1] + det)
+    lambda_m = 1.0 / 2.0 * (gamma_S[0, 0] + gamma_S[1, 1] - det)
+    # compute projectors
+    identity = np.identity(2, np.complex_)
+    c = 1.0 / det
+    e_p = +c * (gamma_S - lambda_m * identity)
+    e_m = -c * (gamma_S - lambda_p * identity)
+    exp = e_m * np.exp(lambda_m) + e_p * np.exp(lambda_p)
+    return exp, lambda_p, lambda_m, e_p, e_m
diff --git a/tests/eko/harmonics/__init__.py b/src/ekore/anomalous_dimensions/polarized/__init__.py
similarity index 100%
rename from tests/eko/harmonics/__init__.py
rename to src/ekore/anomalous_dimensions/polarized/__init__.py
diff --git a/src/ekore/anomalous_dimensions/polarized/space_like/__init__.py b/src/ekore/anomalous_dimensions/polarized/space_like/__init__.py
new file mode 100644
index 000000000..1b64ab235
--- /dev/null
+++ b/src/ekore/anomalous_dimensions/polarized/space_like/__init__.py
@@ -0,0 +1,21 @@
+r"""The polarized, space-like Altarelli-Parisi splitting kernels.
+
+Normalization is given by
+
+.. math::
+    \mathbf{P}(x) = \sum\limits_{j=0} a_s^{j+1} \mathbf P^{(j)}(x)
+
+with :math:`a_s = \frac{\alpha_S(\mu^2)}{4\pi}`.
+"""
+
+import numba as nb
+
+
+@nb.njit(cache=True)
+def gamma_ns(_order, _mode, _n, _nf):
+    raise NotImplementedError("Polarised, space-like is not yet implemented")
+
+
+@nb.njit(cache=True)
+def gamma_singlet(_order, _n, _nf):
+    raise NotImplementedError("Polarised, space-like is not yet implemented")
diff --git a/tests/eko/matching_conditions/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/__init__.py
similarity index 100%
rename from tests/eko/matching_conditions/__init__.py
rename to src/ekore/anomalous_dimensions/unpolarized/__init__.py
diff --git a/src/eko/anomalous_dimensions/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
similarity index 56%
rename from src/eko/anomalous_dimensions/__init__.py
rename to src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
index c3a3f8c78..dc4ad9dbf 100644
--- a/src/eko/anomalous_dimensions/__init__.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
@@ -1,5 +1,4 @@
-r"""
-This module contains the Altarelli-Parisi splitting kernels.
+r"""The unpolarized, space-like Altarelli-Parisi splitting kernels.
 
 Normalization is given by
 
@@ -19,57 +18,8 @@
 import numba as nb
 import numpy as np
 
-from .. import basis_rotation as br
-from .. import harmonics
-from . import aem1, aem2, as1, as1aem1, as2, as3, as4
-
-
-@nb.njit(cache=True)
-def exp_singlet(gamma_S):
-    r"""
-    Computes the exponential and the eigensystem of the singlet anomalous dimension matrix
-
-    Parameters
-    ----------
-        gamma_S : numpy.ndarray
-            singlet anomalous dimension matrix
-
-    Returns
-    -------
-        exp : numpy.ndarray
-            exponential of the singlet anomalous dimension matrix :math:`\gamma_{S}(N)`
-        lambda_p : complex
-            positive eigenvalue of the singlet anomalous dimension matrix
-            :math:`\gamma_{S}(N)`
-        lambda_m : complex
-            negative eigenvalue of the singlet anomalous dimension matrix
-            :math:`\gamma_{S}(N)`
-        e_p : numpy.ndarray
-            projector for the positive eigenvalue of the singlet anomalous
-            dimension matrix :math:`\gamma_{S}(N)`
-        e_m : numpy.ndarray
-            projector for the negative eigenvalue of the singlet anomalous
-            dimension matrix :math:`\gamma_{S}(N)`
-
-    See Also
-    --------
-        eko.anomalous_dimensions.as1.gamma_singlet : :math:`\gamma_{S}^{(0)}(N)`
-        eko.anomalous_dimensions.as2.gamma_singlet : :math:`\gamma_{S}^{(1)}(N)`
-        eko.anomalous_dimensions.as3.gamma_singlet : :math:`\gamma_{S}^{(2)}(N)`
-    """
-    # compute eigenvalues
-    det = np.sqrt(
-        np.power(gamma_S[0, 0] - gamma_S[1, 1], 2) + 4.0 * gamma_S[0, 1] * gamma_S[1, 0]
-    )
-    lambda_p = 1.0 / 2.0 * (gamma_S[0, 0] + gamma_S[1, 1] + det)
-    lambda_m = 1.0 / 2.0 * (gamma_S[0, 0] + gamma_S[1, 1] - det)
-    # compute projectors
-    identity = np.identity(2, np.complex_)
-    c = 1.0 / det
-    e_p = +c * (gamma_S - lambda_m * identity)
-    e_m = -c * (gamma_S - lambda_p * identity)
-    exp = e_m * np.exp(lambda_m) + e_p * np.exp(lambda_p)
-    return exp, lambda_p, lambda_m, e_p, e_m
+from .... import harmonics
+from . import aem1, aem2, as1, as2, as3, as4
 
 
 @nb.njit(cache=True)
@@ -94,15 +44,15 @@ def gamma_ns(order, mode, n, nf):
 
     See Also
     --------
-    eko.anomalous_dimensions.as1.gamma_ns : :math:`\gamma_{ns}^{(0)}(N)`
-    eko.anomalous_dimensions.as2.gamma_nsp : :math:`\gamma_{ns,+}^{(1)}(N)`
-    eko.anomalous_dimensions.as2.gamma_nsm : :math:`\gamma_{ns,-}^{(1)}(N)`
-    eko.anomalous_dimensions.as3.gamma_nsp : :math:`\gamma_{ns,+}^{(2)}(N)`
-    eko.anomalous_dimensions.as3.gamma_nsm : :math:`\gamma_{ns,-}^{(2)}(N)`
-    eko.anomalous_dimensions.as3.gamma_nsv : :math:`\gamma_{ns,v}^{(2)}(N)`
-    eko.anomalous_dimensions.as4.gamma_nsp : :math:`\gamma_{ns,+}^{(3)}(N)`
-    eko.anomalous_dimensions.as4.gamma_nsm : :math:`\gamma_{ns,-}^{(3)}(N)`
-    eko.anomalous_dimensions.as4.gamma_nsv : :math:`\gamma_{ns,v}^{(3)}(N)`
+    ekore.anomalous_dimensions.unpolarized.space_like.as1.gamma_ns : :math:`\gamma_{ns}^{(0)}(N)`
+    ekore.anomalous_dimensions.unpolarized.space_like.as2.gamma_nsp : :math:`\gamma_{ns,+}^{(1)}(N)`
+    ekore.anomalous_dimensions.unpolarized.space_like.as2.gamma_nsm : :math:`\gamma_{ns,-}^{(1)}(N)`
+    ekore.anomalous_dimensions.unpolarized.space_like.as3.gamma_nsp : :math:`\gamma_{ns,+}^{(2)}(N)`
+    ekore.anomalous_dimensions.unpolarized.space_like.as3.gamma_nsm : :math:`\gamma_{ns,-}^{(2)}(N)`
+    ekore.anomalous_dimensions.unpolarized.space_like.as3.gamma_nsv : :math:`\gamma_{ns,v}^{(2)}(N)`
+    ekore.anomalous_dimensions.unpolarized.space_like.as4.gamma_nsp : :math:`\gamma_{ns,+}^{(3)}(N)`
+    ekore.anomalous_dimensions.unpolarized.space_like.as4.gamma_nsm : :math:`\gamma_{ns,-}^{(3)}(N)`
+    ekore.anomalous_dimensions.unpolarized.space_like.as4.gamma_nsv : :math:`\gamma_{ns,v}^{(3)}(N)`
 
     """
     # cache the s-es
@@ -172,10 +122,10 @@ def gamma_singlet(order, n, nf):
 
     See Also
     --------
-    eko.anomalous_dimensions.as1.gamma_singlet : :math:`\gamma_{S}^{(0)}(N)`
-    eko.anomalous_dimensions.as2.gamma_singlet : :math:`\gamma_{S}^{(1)}(N)`
-    eko.anomalous_dimensions.as3.gamma_singlet : :math:`\gamma_{S}^{(2)}(N)`
-    eko.anomalous_dimensions.as4.gamma_singlet : :math:`\gamma_{S}^{(3)}(N)`
+    ekore.anomalous_dimensions.unpolarized.space_like.as1.gamma_singlet : :math:`\gamma_{S}^{(0)}(N)`
+    ekore.anomalous_dimensions.unpolarized.space_like.as2.gamma_singlet : :math:`\gamma_{S}^{(1)}(N)`
+    ekore.anomalous_dimensions.unpolarized.space_like.as3.gamma_singlet : :math:`\gamma_{S}^{(2)}(N)`
+    ekore.anomalous_dimensions.unpolarized.space_like.as4.gamma_singlet : :math:`\gamma_{S}^{(3)}(N)`
 
     """
     # cache the s-es
diff --git a/src/eko/anomalous_dimensions/aem1.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/aem1.py
similarity index 98%
rename from src/eko/anomalous_dimensions/aem1.py
rename to src/ekore/anomalous_dimensions/unpolarized/space_like/aem1.py
index f272439cb..d26683969 100644
--- a/src/eko/anomalous_dimensions/aem1.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/aem1.py
@@ -4,7 +4,7 @@
 
 import numba as nb
 
-from .. import constants
+from eko import constants
 from . import as1
 
 
diff --git a/src/eko/anomalous_dimensions/aem2.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/aem2.py
similarity index 99%
rename from src/eko/anomalous_dimensions/aem2.py
rename to src/ekore/anomalous_dimensions/unpolarized/space_like/aem2.py
index d2a4940d6..83f40e456 100644
--- a/src/eko/anomalous_dimensions/aem2.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/aem2.py
@@ -4,7 +4,7 @@
 
 import numba as nb
 
-from .. import constants
+from eko import constants
 from . import as1aem1
 
 
diff --git a/src/eko/anomalous_dimensions/as1.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as1.py
similarity index 99%
rename from src/eko/anomalous_dimensions/as1.py
rename to src/ekore/anomalous_dimensions/unpolarized/space_like/as1.py
index b4e5fd133..28afb3a00 100644
--- a/src/eko/anomalous_dimensions/as1.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as1.py
@@ -3,7 +3,7 @@
 import numba as nb
 import numpy as np
 
-from .. import constants
+from eko import constants
 
 
 @nb.njit(cache=True)
diff --git a/src/eko/anomalous_dimensions/as1aem1.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as1aem1.py
similarity index 98%
rename from src/eko/anomalous_dimensions/as1aem1.py
rename to src/ekore/anomalous_dimensions/unpolarized/space_like/as1aem1.py
index 85eb25ffa..e16133a05 100644
--- a/src/eko/anomalous_dimensions/as1aem1.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as1aem1.py
@@ -4,8 +4,9 @@
 
 import numba as nb
 
-from .. import constants, harmonics
-from ..harmonics.constants import zeta2, zeta3
+from eko import constants
+from .... import harmonics
+from ....harmonics.constants import zeta2, zeta3
 
 
 @nb.njit(cache=True)
diff --git a/src/eko/anomalous_dimensions/as2.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as2.py
similarity index 99%
rename from src/eko/anomalous_dimensions/as2.py
rename to src/ekore/anomalous_dimensions/unpolarized/space_like/as2.py
index 03945a495..8691a54b2 100644
--- a/src/eko/anomalous_dimensions/as2.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as2.py
@@ -9,8 +9,9 @@
 import numba as nb
 import numpy as np
 
-from .. import constants, harmonics
-from ..harmonics.constants import log2, zeta2, zeta3
+from eko import constants
+from .... import harmonics
+from ....harmonics.constants import log2, zeta2, zeta3
 
 
 @nb.njit(cache=True)
diff --git a/src/eko/anomalous_dimensions/as3.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as3.py
similarity index 99%
rename from src/eko/anomalous_dimensions/as3.py
rename to src/ekore/anomalous_dimensions/unpolarized/space_like/as3.py
index d1b8171a1..a2387bd08 100644
--- a/src/eko/anomalous_dimensions/as3.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as3.py
@@ -8,7 +8,7 @@
 import numba as nb
 import numpy as np
 
-from ..harmonics.constants import zeta2, zeta3
+from ....harmonics.constants import zeta2, zeta3
 
 
 @nb.njit(cache=True)
diff --git a/src/eko/anomalous_dimensions/as4/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py
similarity index 100%
rename from src/eko/anomalous_dimensions/as4/__init__.py
rename to src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py
diff --git a/src/eko/anomalous_dimensions/as4/ggg.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/ggg.py
similarity index 99%
rename from src/eko/anomalous_dimensions/as4/ggg.py
rename to src/ekore/anomalous_dimensions/unpolarized/space_like/as4/ggg.py
index 8570bb3c6..c627b474f 100644
--- a/src/eko/anomalous_dimensions/as4/ggg.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/ggg.py
@@ -3,7 +3,7 @@
 import numba as nb
 import numpy as np
 
-from ...harmonics.log_functions import lm11
+from .....harmonics.log_functions import lm11
 
 
 @nb.njit(cache=True)
diff --git a/src/eko/anomalous_dimensions/as4/ggq.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/ggq.py
similarity index 98%
rename from src/eko/anomalous_dimensions/as4/ggq.py
rename to src/ekore/anomalous_dimensions/unpolarized/space_like/as4/ggq.py
index eb7d96bc1..4752328dc 100644
--- a/src/eko/anomalous_dimensions/as4/ggq.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/ggq.py
@@ -3,7 +3,7 @@
 import numba as nb
 import numpy as np
 
-from ...harmonics.log_functions import lm13, lm13m1, lm14, lm15
+from .....harmonics.log_functions import lm13, lm13m1, lm14, lm15
 
 
 @nb.njit(cache=True)
diff --git a/src/eko/anomalous_dimensions/as4/gnsm.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsm.py
similarity index 97%
rename from src/eko/anomalous_dimensions/as4/gnsm.py
rename to src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsm.py
index 691aede24..336d5b568 100644
--- a/src/eko/anomalous_dimensions/as4/gnsm.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsm.py
@@ -3,9 +3,9 @@
 """
 import numba as nb
 
-from ...constants import CF
-from ...harmonics.constants import zeta3
-from ...harmonics.log_functions import lm11m1, lm12m1, lm13m1
+from eko.constants import CF
+from .....harmonics.constants import zeta3
+from .....harmonics.log_functions import lm11m1, lm12m1, lm13m1
 
 
 @nb.njit(cache=True)
diff --git a/src/eko/anomalous_dimensions/as4/gnsp.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsp.py
similarity index 98%
rename from src/eko/anomalous_dimensions/as4/gnsp.py
rename to src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsp.py
index 1f96efc16..9f5be5c4d 100644
--- a/src/eko/anomalous_dimensions/as4/gnsp.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsp.py
@@ -3,7 +3,7 @@
 """
 import numba as nb
 
-from ...harmonics.log_functions import lm11m1, lm12m1, lm13m1
+from .....harmonics.log_functions import lm11m1, lm12m1, lm13m1
 from .gnsm import gamma_ns_nf3
 
 
diff --git a/src/eko/anomalous_dimensions/as4/gnsv.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsv.py
similarity index 100%
rename from src/eko/anomalous_dimensions/as4/gnsv.py
rename to src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gnsv.py
diff --git a/src/eko/anomalous_dimensions/as4/gps.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gps.py
similarity index 100%
rename from src/eko/anomalous_dimensions/as4/gps.py
rename to src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gps.py
diff --git a/src/eko/anomalous_dimensions/as4/gqg.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gqg.py
similarity index 99%
rename from src/eko/anomalous_dimensions/as4/gqg.py
rename to src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gqg.py
index f23cfef01..7f146eebc 100644
--- a/src/eko/anomalous_dimensions/as4/gqg.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/gqg.py
@@ -3,7 +3,7 @@
 import numba as nb
 import numpy as np
 
-from ...harmonics.log_functions import lm13, lm13m1, lm14, lm15
+from .....harmonics.log_functions import lm13, lm13m1, lm14, lm15
 
 
 @nb.njit(cache=True)
diff --git a/src/ekore/anomalous_dimensions/unpolarized/time_like/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/time_like/__init__.py
new file mode 100644
index 000000000..7d353c93e
--- /dev/null
+++ b/src/ekore/anomalous_dimensions/unpolarized/time_like/__init__.py
@@ -0,0 +1,21 @@
+r"""The unpolarized, time-like Altarelli-Parisi splitting kernels.
+
+Normalization is given by
+
+.. math::
+    \mathbf{P}(x) = \sum\limits_{j=0} a_s^{j+1} \mathbf P^{(j)}(x)
+
+with :math:`a_s = \frac{\alpha_S(\mu^2)}{4\pi}`.
+"""
+
+import numba as nb
+
+
+@nb.njit(cache=True)
+def gamma_ns(_order, _mode, _n, _nf):
+    raise NotImplementedError("Polarised is not yet implemented")
+
+
+@nb.njit(cache=True)
+def gamma_singlet(_order, _n, _nf):
+    raise NotImplementedError("Polarised is not yet implemented")
diff --git a/src/eko/harmonics/__init__.py b/src/ekore/harmonics/__init__.py
similarity index 100%
rename from src/eko/harmonics/__init__.py
rename to src/ekore/harmonics/__init__.py
diff --git a/src/eko/harmonics/constants.py b/src/ekore/harmonics/constants.py
similarity index 100%
rename from src/eko/harmonics/constants.py
rename to src/ekore/harmonics/constants.py
diff --git a/src/eko/harmonics/f_functions/__init__.py b/src/ekore/harmonics/f_functions/__init__.py
similarity index 100%
rename from src/eko/harmonics/f_functions/__init__.py
rename to src/ekore/harmonics/f_functions/__init__.py
diff --git a/src/eko/harmonics/f_functions/f11.py b/src/ekore/harmonics/f_functions/f11.py
similarity index 100%
rename from src/eko/harmonics/f_functions/f11.py
rename to src/ekore/harmonics/f_functions/f11.py
diff --git a/src/eko/harmonics/f_functions/f13.py b/src/ekore/harmonics/f_functions/f13.py
similarity index 100%
rename from src/eko/harmonics/f_functions/f13.py
rename to src/ekore/harmonics/f_functions/f13.py
diff --git a/src/eko/harmonics/f_functions/f14_f12.py b/src/ekore/harmonics/f_functions/f14_f12.py
similarity index 100%
rename from src/eko/harmonics/f_functions/f14_f12.py
rename to src/ekore/harmonics/f_functions/f14_f12.py
diff --git a/src/eko/harmonics/f_functions/f16.py b/src/ekore/harmonics/f_functions/f16.py
similarity index 100%
rename from src/eko/harmonics/f_functions/f16.py
rename to src/ekore/harmonics/f_functions/f16.py
diff --git a/src/eko/harmonics/f_functions/f17.py b/src/ekore/harmonics/f_functions/f17.py
similarity index 100%
rename from src/eko/harmonics/f_functions/f17.py
rename to src/ekore/harmonics/f_functions/f17.py
diff --git a/src/eko/harmonics/f_functions/f18.py b/src/ekore/harmonics/f_functions/f18.py
similarity index 100%
rename from src/eko/harmonics/f_functions/f18.py
rename to src/ekore/harmonics/f_functions/f18.py
diff --git a/src/eko/harmonics/f_functions/f19.py b/src/ekore/harmonics/f_functions/f19.py
similarity index 100%
rename from src/eko/harmonics/f_functions/f19.py
rename to src/ekore/harmonics/f_functions/f19.py
diff --git a/src/eko/harmonics/f_functions/f20.py b/src/ekore/harmonics/f_functions/f20.py
similarity index 100%
rename from src/eko/harmonics/f_functions/f20.py
rename to src/ekore/harmonics/f_functions/f20.py
diff --git a/src/eko/harmonics/f_functions/f21.py b/src/ekore/harmonics/f_functions/f21.py
similarity index 100%
rename from src/eko/harmonics/f_functions/f21.py
rename to src/ekore/harmonics/f_functions/f21.py
diff --git a/src/eko/harmonics/f_functions/f9.py b/src/ekore/harmonics/f_functions/f9.py
similarity index 100%
rename from src/eko/harmonics/f_functions/f9.py
rename to src/ekore/harmonics/f_functions/f9.py
diff --git a/src/eko/harmonics/g_functions.py b/src/ekore/harmonics/g_functions.py
similarity index 100%
rename from src/eko/harmonics/g_functions.py
rename to src/ekore/harmonics/g_functions.py
diff --git a/src/eko/harmonics/log_functions.py b/src/ekore/harmonics/log_functions.py
similarity index 100%
rename from src/eko/harmonics/log_functions.py
rename to src/ekore/harmonics/log_functions.py
diff --git a/src/eko/harmonics/polygamma.py b/src/ekore/harmonics/polygamma.py
similarity index 100%
rename from src/eko/harmonics/polygamma.py
rename to src/ekore/harmonics/polygamma.py
diff --git a/src/eko/harmonics/w1.py b/src/ekore/harmonics/w1.py
similarity index 95%
rename from src/eko/harmonics/w1.py
rename to src/ekore/harmonics/w1.py
index 5d4e5b27c..10f7cc3d4 100644
--- a/src/eko/harmonics/w1.py
+++ b/src/ekore/harmonics/w1.py
@@ -28,7 +28,7 @@ def S1(N):
 
     See Also
     --------
-        eko.harmonics.polygamma.cern_polygamma : :math:`\psi_k(N)`
+        ekore.harmonics.polygamma.cern_polygamma : :math:`\psi_k(N)`
     """
     return cern_polygamma(N + 1.0, 0) + np.euler_gamma
 
diff --git a/src/eko/harmonics/w2.py b/src/ekore/harmonics/w2.py
similarity index 95%
rename from src/eko/harmonics/w2.py
rename to src/ekore/harmonics/w2.py
index 9e978f0d4..69c0029eb 100644
--- a/src/eko/harmonics/w2.py
+++ b/src/ekore/harmonics/w2.py
@@ -28,7 +28,7 @@ def S2(N):
 
     See Also
     --------
-        eko.harmonics.polygamma.cern_polygamma : :math:`\psi_k(N)`
+        ekore.harmonics.polygamma.cern_polygamma : :math:`\psi_k(N)`
     """
     return -cern_polygamma(N + 1.0, 1) + zeta2
 
diff --git a/src/eko/harmonics/w3.py b/src/ekore/harmonics/w3.py
similarity index 93%
rename from src/eko/harmonics/w3.py
rename to src/ekore/harmonics/w3.py
index 49e3d5a7e..dcd52428c 100644
--- a/src/eko/harmonics/w3.py
+++ b/src/ekore/harmonics/w3.py
@@ -29,7 +29,7 @@ def S3(N):
 
     See Also
     --------
-    eko.harmonics.polygamma.cern_polygamma : :math:`\psi_k(N)`
+    ekore.harmonics.polygamma.cern_polygamma : :math:`\psi_k(N)`
 
     """
     return 0.5 * cern_polygamma(N + 1.0, 2) + zeta3
@@ -59,7 +59,7 @@ def Sm3(N, hS3, is_singlet=None):
 
     See Also
     --------
-    eko.harmonics.w3.S3 : :math:`S_3(N)`
+    ekore.harmonics.w3.S3 : :math:`S_3(N)`
 
     """
     if is_singlet is None:
@@ -97,7 +97,7 @@ def S21(N, S1, S2):
 
     See Also
     --------
-    eko.harmonics.g_functions.mellin_g18 : :math:`g_18(N)`
+    ekore.harmonics.g_functions.mellin_g18 : :math:`g_18(N)`
 
     """
     return -gf.mellin_g18(N, S1, S2) + 2 * zeta3
@@ -129,7 +129,7 @@ def Sm21(N, S1, Sm1, is_singlet=None):
 
     See Also
     --------
-    eko.harmonics.g_functions : :math:`g_3(N)`
+    ekore.harmonics.g_functions : :math:`g_3(N)`
 
     """
     # Note mellin g3 was integrated following x^(N-1) convention.
@@ -170,7 +170,7 @@ def S2m1(N, S2, Sm1, Sm2, is_singlet=None):
 
     See Also
     --------
-    eko.harmonics.g_functions.mellin_g4 : :math:`g_4(N)`
+    ekore.harmonics.g_functions.mellin_g4 : :math:`g_4(N)`
 
     """
     eta = symmetry_factor(N, is_singlet)
@@ -208,7 +208,7 @@ def Sm2m1(N, S1, S2, Sm2):
 
     See Also
     --------
-    eko.harmonics.g_functions.mellin_g19 : :math:`g_19(N)`
+    ekore.harmonics.g_functions.mellin_g19 : :math:`g_19(N)`
 
     """
     return -gf.mellin_g19(N, S1) + log2 * (S2 - Sm2) - 5 / 8 * zeta3
diff --git a/src/eko/harmonics/w4.py b/src/ekore/harmonics/w4.py
similarity index 94%
rename from src/eko/harmonics/w4.py
rename to src/ekore/harmonics/w4.py
index 270e79311..855602ec7 100644
--- a/src/eko/harmonics/w4.py
+++ b/src/ekore/harmonics/w4.py
@@ -29,7 +29,7 @@ def S4(N):
 
     See Also
     --------
-    eko.harmonics.polygamma.cern_polygamma : :math:`\psi_k(N)`
+    ekore.harmonics.polygamma.cern_polygamma : :math:`\psi_k(N)`
 
     """
     return zeta4 - 1.0 / 6.0 * cern_polygamma(N + 1.0, 3)
@@ -102,7 +102,7 @@ def Sm31(N, S1, Sm1, Sm2, is_singlet=None):
 
     See Also
     --------
-    eko.harmonics.g_functions.mellin_g6 : :math:`g_6(N)`
+    ekore.harmonics.g_functions.mellin_g6 : :math:`g_6(N)`
 
     """
     eta = symmetry_factor(N, is_singlet)
@@ -148,7 +148,7 @@ def Sm22(N, S1, S2, Sm2, Sm31, is_singlet=None):
 
     See Also
     --------
-    eko.harmonics.g_functions.mellin_g5 : :math:`g_5(N)`
+    ekore.harmonics.g_functions.mellin_g5 : :math:`g_5(N)`
 
     """
     eta = symmetry_factor(N, is_singlet)
@@ -185,7 +185,7 @@ def Sm211(N, S1, S2, Sm1, is_singlet=None):
 
     See Also
     --------
-    eko.harmonics.g_functions.mellin_g8 : :math:`g_8(N)`
+    ekore.harmonics.g_functions.mellin_g8 : :math:`g_8(N)`
 
     """
     eta = symmetry_factor(N, is_singlet)
@@ -225,7 +225,7 @@ def S211(N, S1, S2, S3):
 
     See Also
     --------
-    eko.harmonics.g_functions.mellin_g21 : :math:`g_21(N)`
+    ekore.harmonics.g_functions.mellin_g21 : :math:`g_21(N)`
 
     """
     return -gf.mellin_g21(N, S1, S2, S3) + 6 / 5 * zeta2**2
@@ -258,7 +258,7 @@ def S31(N, S1, S2, S3, S4):
 
     See Also
     --------
-    eko.harmonics.g_functions.mellin_g22 : :math:`g_22(N)`
+    ekore.harmonics.g_functions.mellin_g22 : :math:`g_22(N)`
 
     """
     return (
diff --git a/src/eko/harmonics/w5.py b/src/ekore/harmonics/w5.py
similarity index 99%
rename from src/eko/harmonics/w5.py
rename to src/ekore/harmonics/w5.py
index d6a3f91d9..b75f9eebf 100644
--- a/src/eko/harmonics/w5.py
+++ b/src/ekore/harmonics/w5.py
@@ -29,7 +29,7 @@ def S5(N):
 
     See Also
     --------
-    eko.harmonics.polygamma.cern_polygamma : :math:`\psi_k(N)`
+    ekore.harmonics.polygamma.cern_polygamma : :math:`\psi_k(N)`
 
     """
     return zeta5 + 1.0 / 24.0 * cern_polygamma(N + 1.0, 4)
diff --git a/src/ekore/operator_matrix_elements/__init__.py b/src/ekore/operator_matrix_elements/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/src/ekore/operator_matrix_elements/polarized/__init__.py b/src/ekore/operator_matrix_elements/polarized/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/src/ekore/operator_matrix_elements/polarized/space_like/__init__.py b/src/ekore/operator_matrix_elements/polarized/space_like/__init__.py
new file mode 100644
index 000000000..5d5f771f7
--- /dev/null
+++ b/src/ekore/operator_matrix_elements/polarized/space_like/__init__.py
@@ -0,0 +1,13 @@
+r"""The polarized, space-like |OME|."""
+
+import numba as nb
+
+
+@nb.njit(cache=True)
+def A_non_singlet(_matching_order, _n, _sx, _nf, _L):
+    raise NotImplementedError("Polarised, space-like is not yet implemented")
+
+
+@nb.njit(cache=True)
+def A_singlet(_matching_order, _n, _sx, _nf, _L, _is_msbar, _sx_ns=None):
+    raise NotImplementedError("Polarised, space-like is not yet implemented")
diff --git a/src/ekore/operator_matrix_elements/unpolarized/__init__.py b/src/ekore/operator_matrix_elements/unpolarized/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/src/ekore/operator_matrix_elements/unpolarized/space_like/__init__.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/__init__.py
new file mode 100644
index 000000000..a1d2f8ec6
--- /dev/null
+++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/__init__.py
@@ -0,0 +1,91 @@
+"""
+This module defines the |OME| for the non-trivial matching conditions in the
+|VFNS| evolution.
+"""
+
+import numba as nb
+import numpy as np
+
+from . import as1, as2, as3
+
+
+@nb.njit(cache=True)
+def A_singlet(matching_order, n, sx, nf, L, is_msbar, sx_ns=None):
+    r"""
+    Computes the tower of the singlet |OME|.
+
+    Parameters
+    ----------
+        matching_order : tuple(int,int)
+            perturbative matching_order
+        n : complex
+            Mellin variable
+        sx : list
+            singlet like harmonic sums cache
+        nf: int
+            number of active flavor below threshold
+        L : float
+            :math:``\ln(\mu_F^2 / m_h^2)``
+        is_msbar: bool
+            add the |MSbar| contribution
+        sx_ns : list
+            non-singlet like harmonic sums cache
+
+    Returns
+    -------
+        A_singlet : numpy.ndarray
+            singlet |OME|
+
+    See Also
+    --------
+        ekore.matching_conditions.nlo.A_singlet_1 : :math:`A^{S,(1)}(N)`
+        ekore.matching_conditions.nlo.A_hh_1 : :math:`A_{HH}^{(1)}(N)`
+        ekore.matching_conditions.nlo.A_gh_1 : :math:`A_{gH}^{(1)}(N)`
+        ekore.matching_conditions.nnlo.A_singlet_2 : :math:`A_{S,(2)}(N)`
+    """
+    A_s = np.zeros((matching_order[0], 3, 3), np.complex_)
+    if matching_order[0] >= 1:
+        A_s[0] = as1.A_singlet(n, sx, L)
+    if matching_order[0] >= 2:
+        A_s[1] = as2.A_singlet(n, sx, L, is_msbar)
+    if matching_order[0] >= 3:
+        A_s[2] = as3.A_singlet(n, sx, sx_ns, nf, L)
+    return A_s
+
+
+@nb.njit(cache=True)
+def A_non_singlet(matching_order, n, sx, nf, L):
+    r"""
+    Computes the tower of the non-singlet |OME|
+
+    Parameters
+    ----------
+        matching_order : tuple(int,int)
+            perturbative matching_order
+        n : complex
+            Mellin variable
+        sx : list
+            harmonic sums cache
+        nf: int
+            number of active flavor below threshold
+        L : float
+            :math:``\ln(\mu_F^2 / m_h^2)``
+
+    Returns
+    -------
+        A_non_singlet : numpy.ndarray
+            non-singlet |OME|
+
+    See Also
+    --------
+        ekore.matching_conditions.nlo.A_hh_1 : :math:`A_{HH}^{(1)}(N)`
+        ekore.matching_conditions.nnlo.A_ns_2 : :math:`A_{qq,H}^{NS,(2)}`
+    """
+    A_ns = np.zeros((matching_order[0], 2, 2), np.complex_)
+    if matching_order[0] >= 1:
+        A_ns[0] = as1.A_ns(n, sx, L)
+    if matching_order[0] >= 2:
+        A_ns[1] = as2.A_ns(n, sx, L)
+    if matching_order[0] >= 3:
+        A_ns[2] = as3.A_ns(n, sx, nf, L)
+    return A_ns
diff --git a/src/eko/matching_conditions/as1.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as1.py
similarity index 99%
rename from src/eko/matching_conditions/as1.py
rename to src/ekore/operator_matrix_elements/unpolarized/space_like/as1.py
index cc74550ca..c373f8091 100644
--- a/src/eko/matching_conditions/as1.py
+++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as1.py
@@ -9,7 +9,7 @@
 import numba as nb
 import numpy as np
 
-from ..constants import CF
+from eko.constants import CF
 
 
 @nb.njit(cache=True)
diff --git a/src/eko/matching_conditions/as2.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as2.py
similarity index 99%
rename from src/eko/matching_conditions/as2.py
rename to src/ekore/operator_matrix_elements/unpolarized/space_like/as2.py
index 7a86d8050..2197fbfd9 100644
--- a/src/eko/matching_conditions/as2.py
+++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as2.py
@@ -11,8 +11,9 @@
 import numba as nb
 import numpy as np
 
-from .. import constants
-from ..harmonics.constants import zeta2, zeta3
+from eko import constants
+
+from ....harmonics.constants import zeta2, zeta3
 from .as1 import A_gg as A_gg_1
 from .as1 import A_hg as A_hg_1
 
diff --git a/src/eko/matching_conditions/as3/__init__.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/__init__.py
similarity index 97%
rename from src/eko/matching_conditions/as3/__init__.py
rename to src/ekore/operator_matrix_elements/unpolarized/space_like/as3/__init__.py
index b414369a5..06feccd38 100644
--- a/src/eko/matching_conditions/as3/__init__.py
+++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/__init__.py
@@ -80,7 +80,7 @@ def A_singlet(n, sx_singlet, sx_non_singlet, nf, L):
         \end{array}\right)
 
     When using the code, please cite the complete list of references
-    available at the top of this module :mod:`eko.matching_conditions.as3`.
+    available at the top of this module :mod:`ekore.matching_conditions.as3`.
 
     Parameters
     ----------
@@ -140,7 +140,7 @@ def A_ns(n, sx_all, nf, L):
         \end{array}\right)
 
     When using the code, please cite the complete list of references available
-    at the top of this module :mod:`eko.matching_conditions.as3`.
+    at the top of this module :mod:`ekore.matching_conditions.as3`.
 
     Parameters
     ----------
diff --git a/src/eko/matching_conditions/as3/aHg.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHg.py
similarity index 99%
rename from src/eko/matching_conditions/as3/aHg.py
rename to src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHg.py
index b756ab584..c676ce558 100644
--- a/src/eko/matching_conditions/as3/aHg.py
+++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHg.py
@@ -11,7 +11,7 @@ def A_Hg(n, sx, nf, L):  # pylint: disable=too-many-locals
     The expression is presented in :cite:`Bierenbaum:2009mv`.
 
     When using the code, please cite the complete list of references
-    available in :mod:`eko.matching_conditions.as3`.
+    available in :mod:`ekore.matching_conditions.as3`.
 
     Parameters
     ----------
@@ -31,7 +31,7 @@ def A_Hg(n, sx, nf, L):  # pylint: disable=too-many-locals
 
     See Also
     --------
-    A_Hgstfac: eko.matching_conditions.as3.aHgstfac.A_Hgstfac
+    A_Hgstfac: ekore.matching_conditions.as3.aHgstfac.A_Hgstfac
         Incomplete part of the |OME|.
 
     """
diff --git a/src/eko/matching_conditions/as3/aHgstfac.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHgstfac.py
similarity index 99%
rename from src/eko/matching_conditions/as3/aHgstfac.py
rename to src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHgstfac.py
index 863c46536..4eed1890a 100644
--- a/src/eko/matching_conditions/as3/aHgstfac.py
+++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHgstfac.py
@@ -10,7 +10,7 @@ def A_Hgstfac(n, sx, nf):
     The expression is presented in :cite:`Blumlein:2017wxd` (eq 3.1).
 
     When using the code, please cite the complete list of references
-    available in :mod:`eko.matching_conditions.as3`.
+    available in :mod:`ekore.matching_conditions.as3`.
 
     The expression contains some unknown parts which are set to 0.
     However we have included a shift to impose the gluon
diff --git a/src/eko/matching_conditions/as3/aHq.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHq.py
similarity index 99%
rename from src/eko/matching_conditions/as3/aHq.py
rename to src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHq.py
index 18a0bb61c..6ea2968ec 100644
--- a/src/eko/matching_conditions/as3/aHq.py
+++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aHq.py
@@ -10,7 +10,7 @@ def A_Hq(n, sx, nf, L):  # pylint: disable=too-many-locals
     and :cite:`Blumlein:2017wxd` (eq 3.1).
 
     When using the code, please cite the complete list of references
-    available in :mod:`eko.matching_conditions.as3`.
+    available in :mod:`ekore.matching_conditions.as3`.
 
     The part proportional to :math:`n_f^0` includes non trivial weight-5
     harmonics and has been parametrized in Mellin space.
diff --git a/src/eko/matching_conditions/as3/agg.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/agg.py
similarity index 99%
rename from src/eko/matching_conditions/as3/agg.py
rename to src/ekore/operator_matrix_elements/unpolarized/space_like/as3/agg.py
index 79debe7c3..b4783d83e 100644
--- a/src/eko/matching_conditions/as3/agg.py
+++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/agg.py
@@ -11,7 +11,7 @@ def A_gg(n, sx, nf, L):  # pylint: disable=too-many-locals
     The expression is presented in :cite:`Bierenbaum:2009mv`.
 
     When using the code, please cite the complete list of references
-    available in :mod:`eko.matching_conditions.as3`.
+    available in :mod:`ekore.matching_conditions.as3`.
 
     Parameters
     ----------
@@ -31,7 +31,7 @@ def A_gg(n, sx, nf, L):  # pylint: disable=too-many-locals
 
     See Also
     --------
-    A_ggTF2: eko.matching_conditions.as3.aggTF2.A_ggTF2
+    A_ggTF2: ekore.matching_conditions.as3.aggTF2.A_ggTF2
         Incomplete part proportional to :math:`T_{F}^2`.
 
     """
diff --git a/src/eko/matching_conditions/as3/aggTF2.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aggTF2.py
similarity index 99%
rename from src/eko/matching_conditions/as3/aggTF2.py
rename to src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aggTF2.py
index 5e337c5f0..299006dc6 100644
--- a/src/eko/matching_conditions/as3/aggTF2.py
+++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aggTF2.py
@@ -10,7 +10,7 @@ def A_ggTF2(n, sx):
     It contains a binomial factor which is given approximated.
 
     When using the code, please cite the complete list of references
-    available in :mod:`eko.matching_conditions.as3`.
+    available in :mod:`ekore.matching_conditions.as3`.
 
     Parameters
     ----------
diff --git a/src/eko/matching_conditions/as3/agq.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/agq.py
similarity index 99%
rename from src/eko/matching_conditions/as3/agq.py
rename to src/ekore/operator_matrix_elements/unpolarized/space_like/as3/agq.py
index 1d874021d..906e08e52 100644
--- a/src/eko/matching_conditions/as3/agq.py
+++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/agq.py
@@ -8,7 +8,7 @@ def A_gq(n, sx, nf, L):  # pylint: disable=too-many-locals
     The expression is presented in :cite:`Ablinger_2014` (eq 6.3).
 
     When using the code, please cite the complete list of references
-    available in :mod:`eko.matching_conditions.as3`.
+    available in :mod:`ekore.matching_conditions.as3`.
 
     Parameters
     ----------
diff --git a/src/eko/matching_conditions/as3/aqg.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqg.py
similarity index 99%
rename from src/eko/matching_conditions/as3/aqg.py
rename to src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqg.py
index 33a059eb9..3c2cc36e4 100644
--- a/src/eko/matching_conditions/as3/aqg.py
+++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqg.py
@@ -8,7 +8,7 @@ def A_qg(n, sx, nf, L):
     The expression is presented in :cite:`Bierenbaum:2009mv`.
 
     When using the code, please cite the complete list of references
-    available in :mod:`eko.matching_conditions.as3`.
+    available in :mod:`ekore.matching_conditions.as3`.
 
     Parameters
     ----------
diff --git a/src/eko/matching_conditions/as3/aqqNS.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqqNS.py
similarity index 99%
rename from src/eko/matching_conditions/as3/aqqNS.py
rename to src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqqNS.py
index fdbf9114e..40be9ed69 100644
--- a/src/eko/matching_conditions/as3/aqqNS.py
+++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqqNS.py
@@ -9,7 +9,7 @@ def A_qqNS(n, sx, nf, L):
     :cite:`Ablinger:2014vwa`. It contains some weight 5 harmonics sums.
 
     When using the code, please cite the complete list of references
-    available in :mod:`eko.matching_conditions.as3`.
+    available in :mod:`ekore.matching_conditions.as3`.
 
     Note the part proportional to nf^0 includes weight = 5
     harmonics and has been parametrized in Mellin space.
diff --git a/src/eko/matching_conditions/as3/aqqPS.py b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqqPS.py
similarity index 99%
rename from src/eko/matching_conditions/as3/aqqPS.py
rename to src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqqPS.py
index a0f20af1c..ee8f136cd 100644
--- a/src/eko/matching_conditions/as3/aqqPS.py
+++ b/src/ekore/operator_matrix_elements/unpolarized/space_like/as3/aqqPS.py
@@ -8,7 +8,7 @@ def A_qqPS(n, sx, nf, L):
     The expression is presented in :cite:`Bierenbaum:2009mv`.
 
     When using the code, please cite the complete list of references
-    available in :mod:`eko.matching_conditions.as3`.
+    available in :mod:`ekore.matching_conditions.as3`.
 
     Parameters
     ----------
diff --git a/src/ekore/operator_matrix_elements/unpolarized/time_like/__init__.py b/src/ekore/operator_matrix_elements/unpolarized/time_like/__init__.py
new file mode 100644
index 000000000..53c68e8d6
--- /dev/null
+++ b/src/ekore/operator_matrix_elements/unpolarized/time_like/__init__.py
@@ -0,0 +1,13 @@
+r"""The polarized, time-like |OME|."""
+
+import numba as nb
+
+
+@nb.njit(cache=True)
+def A_non_singlet(_matching_order, _n, _sx, _nf, _L):
+    raise NotImplementedError("Time-like is not yet implemented")
+
+
+@nb.njit(cache=True)
+def A_singlet(_matching_order, _n, _sx, _nf, _L, _is_msbar, _sx_ns=None):
+    raise NotImplementedError("Time-like is not yet implemented")
diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py
index d40fbeb39..c33aa13a4 100644
--- a/tests/eko/evolution_operator/test_init.py
+++ b/tests/eko/evolution_operator/test_init.py
@@ -4,7 +4,7 @@
 import scipy.integrate
 
 import eko.runner.legacy
-from eko import anomalous_dimensions as ad
+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
@@ -43,6 +43,8 @@ def test_quad_ker(monkeypatch):
             ev_op_max_order=(0, 0),
             sv_mode=1,
             is_threshold=False,
+            is_polarized=False,
+            is_time_like=False,
         )
         np.testing.assert_allclose(res_ns, 0.0)
         res_s = quad_ker(
@@ -62,6 +64,8 @@ def test_quad_ker(monkeypatch):
             ev_op_max_order=(0, 0),
             sv_mode=1,
             is_threshold=False,
+            is_polarized=False,
+            is_time_like=False,
         )
         np.testing.assert_allclose(res_s, 1.0)
         res_s = quad_ker(
@@ -81,6 +85,8 @@ def test_quad_ker(monkeypatch):
             ev_op_max_order=(0, 0),
             sv_mode=1,
             is_threshold=False,
+            is_polarized=False,
+            is_time_like=False,
         )
         np.testing.assert_allclose(res_s, 0.0)
     for label in [(br.non_singlet_pids_map["ns+"], 0), (100, 100)]:
@@ -102,6 +108,8 @@ def test_quad_ker(monkeypatch):
                 ev_op_max_order=(1, 0),
                 sv_mode=sv,
                 is_threshold=False,
+                is_polarized=False,
+                is_time_like=False,
             )
             np.testing.assert_allclose(res_sv, 1.0)
 
@@ -123,6 +131,8 @@ def test_quad_ker(monkeypatch):
         ev_op_max_order=(0, 0),
         sv_mode=1,
         is_threshold=False,
+        is_polarized=False,
+        is_time_like=False,
     )
     np.testing.assert_allclose(res_ns, 0.0)
 
@@ -336,6 +346,8 @@ def quad_ker_pegasus(
                     10,
                     0,
                     False,
+                    False,
+                    False,
                 ),
                 epsabs=1e-12,
                 epsrel=1e-5,
diff --git a/tests/eko/matching_conditions/test_init.py b/tests/eko/evolution_operator/test_matchin_condition.py
similarity index 98%
rename from tests/eko/matching_conditions/test_init.py
rename to tests/eko/evolution_operator/test_matchin_condition.py
index ab93c3487..cf6ddd0f5 100644
--- a/tests/eko/matching_conditions/test_init.py
+++ b/tests/eko/evolution_operator/test_matchin_condition.py
@@ -3,7 +3,7 @@
 
 from eko import basis_rotation as br
 from eko import member
-from eko.matching_conditions import MatchingCondition
+from eko.evolution_operator.matching_condition import MatchingCondition
 
 
 def mkOM(shape):
diff --git a/tests/eko/matching_conditions/test_ome.py b/tests/eko/evolution_operator/test_ome.py
similarity index 95%
rename from tests/eko/matching_conditions/test_ome.py
rename to tests/eko/evolution_operator/test_ome.py
index 2df26a03f..0d33d071d 100644
--- a/tests/eko/matching_conditions/test_ome.py
+++ b/tests/eko/evolution_operator/test_ome.py
@@ -1,21 +1,23 @@
-# Test eko.matching_conditions.OperatorMatrixElement
+# Test ekore.matching_conditions.OperatorMatrixElement
 import pathlib
 
 import numpy as np
 
 from eko import basis_rotation as br
 from eko import interpolation, mellin
-from eko.harmonics import compute_cache
-from eko.io.runcards import OperatorCard, TheoryCard
-from eko.io.types import InversionMethod
-from eko.matching_conditions.operator_matrix_element import (
-    A_non_singlet,
-    A_singlet,
+from eko.evolution_operator.operator_matrix_element import (
     OperatorMatrixElement,
     build_ome,
     quad_ker,
 )
+from eko.io.runcards import OperatorCard, TheoryCard
+from eko.io.types import InversionMethod
 from eko.runner import legacy
+from ekore.harmonics import compute_cache
+from ekore.operator_matrix_elements.unpolarized.space_like import (
+    A_non_singlet,
+    A_singlet,
+)
 
 max_weight_dict = {1: 2, 2: 3, 3: 5}
 
@@ -88,12 +90,12 @@ def test_quad_ker(monkeypatch):
     monkeypatch.setattr(interpolation, "evaluate_Nx", lambda *args: 1)
     zeros = np.zeros((2, 2))
     monkeypatch.setattr(
-        "eko.matching_conditions.operator_matrix_element.A_non_singlet",
+        "ekore.operator_matrix_elements.unpolarized.space_like.A_non_singlet",
         lambda *args: np.array([zeros, zeros, zeros]),
     )
     zeros = np.zeros((3, 3))
     monkeypatch.setattr(
-        "eko.matching_conditions.operator_matrix_element.A_singlet",
+        "ekore.operator_matrix_elements.unpolarized.space_like.A_singlet",
         lambda *args: np.array([zeros, zeros, zeros]),
     )
     for is_log in [True, False]:
@@ -110,6 +112,8 @@ def test_quad_ker(monkeypatch):
             nf=3,
             L=0.0,
             is_msbar=False,
+            is_polarized=False,
+            is_time_like=False,
         )
         np.testing.assert_allclose(res_ns, 1.0)
         res_s = quad_ker(
@@ -125,6 +129,8 @@ def test_quad_ker(monkeypatch):
             nf=3,
             L=0.0,
             is_msbar=False,
+            is_polarized=False,
+            is_time_like=False,
         )
         np.testing.assert_allclose(res_s, 1.0)
         res_s = quad_ker(
@@ -140,6 +146,8 @@ def test_quad_ker(monkeypatch):
             nf=3,
             L=0.0,
             is_msbar=False,
+            is_polarized=False,
+            is_time_like=False,
         )
         np.testing.assert_allclose(res_s, 0.0)
 
@@ -159,6 +167,8 @@ def test_quad_ker(monkeypatch):
             nf=3,
             L=0.0,
             is_msbar=False,
+            is_polarized=False,
+            is_time_like=False,
         )
         if label[-1] == label[-2]:
             np.testing.assert_allclose(res_ns, 1.0)
@@ -192,6 +202,8 @@ def test_quad_ker(monkeypatch):
             nf=3,
             L=0.0,
             is_msbar=False,
+            is_polarized=False,
+            is_time_like=False,
         )
         if label[-1] == label[-2]:
             np.testing.assert_allclose(res_ns, 1.0)
@@ -212,6 +224,8 @@ def test_quad_ker(monkeypatch):
         nf=3,
         L=0.0,
         is_msbar=False,
+        is_polarized=False,
+        is_time_like=False,
     )
     np.testing.assert_allclose(res_ns, 0.0)
 
diff --git a/tests/eko/evolution_operator/test_physical.py b/tests/eko/evolution_operator/test_physical.py
index a02591d44..a194f0502 100644
--- a/tests/eko/evolution_operator/test_physical.py
+++ b/tests/eko/evolution_operator/test_physical.py
@@ -3,8 +3,8 @@
 
 from eko import basis_rotation as br
 from eko import member
+from eko.evolution_operator.matching_condition import MatchingCondition
 from eko.evolution_operator.physical import PhysicalOperator
-from eko.matching_conditions import MatchingCondition
 
 
 def mkOM(shape):
diff --git a/tests/eko/kernels/test_ns.py b/tests/eko/kernels/test_ns.py
index 17d4964b8..9e89492aa 100644
--- a/tests/eko/kernels/test_ns.py
+++ b/tests/eko/kernels/test_ns.py
@@ -3,7 +3,7 @@
 import numpy as np
 import pytest
 
-from eko import anomalous_dimensions as ad
+import ekore.anomalous_dimensions.unpolarized.space_like as ad
 from eko import beta
 from eko.kernels import non_singlet as ns
 
diff --git a/tests/eko/kernels/test_s.py b/tests/eko/kernels/test_s.py
index ec4d93293..618699e84 100644
--- a/tests/eko/kernels/test_s.py
+++ b/tests/eko/kernels/test_s.py
@@ -3,7 +3,7 @@
 import numpy as np
 import pytest
 
-from eko import anomalous_dimensions as ad
+from ekore import anomalous_dimensions as ad
 from eko.kernels import singlet as s
 
 methods = [
diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py
index c1ecc33fd..2a95edcb3 100644
--- a/tests/eko/scale_variations/test_expanded.py
+++ b/tests/eko/scale_variations/test_expanded.py
@@ -1,7 +1,7 @@
 import numpy as np
 
 from eko import basis_rotation as br
-from eko.anomalous_dimensions import gamma_ns, gamma_singlet
+from ekore.anomalous_dimensions.unpolarized.space_like import gamma_ns, gamma_singlet
 from eko.beta import beta_qcd_as2, beta_qcd_as3
 from eko.kernels import non_singlet, singlet
 from eko.scale_variations import Modes, expanded, exponentiated
@@ -86,7 +86,7 @@ def scheme_diff(g, k, pto, is_singlet):
                 - 1 / 2 * a0**3 * b1 * g[0] * k**2
                 + a1**3 * b1 * g[0] * k**2
                 - a0**3 * b0 * g[1] * k**2
-                + 2 * a1**3 *b0 * g[1] * k**2
+                + 2 * a1**3 * b0 * g[1] * k**2
                 + a0**3 * g0g1 * k**2
                 - a0**2 * a1 * g0g1 * k**2
                 - a0 * a1**2 * g0g1 * k**2
@@ -104,9 +104,8 @@ def scheme_diff(g, k, pto, is_singlet):
             )
         return diff
 
-
     for L in [np.log(0.5), np.log(2)]:
-        for order in [(2, 0), (3, 0), (4,0)]:
+        for order in [(2, 0), (3, 0), (4, 0)]:
             # Non singlet kernels
             gns = gamma_ns(order, br.non_singlet_pids_map["ns+"], n, nf)
             ker = non_singlet.dispatcher(
diff --git a/tests/eko/test_beta.py b/tests/eko/test_beta.py
index 4e2eb814e..fb089bb62 100644
--- a/tests/eko/test_beta.py
+++ b/tests/eko/test_beta.py
@@ -6,7 +6,7 @@
 import pytest
 
 from eko import beta
-from eko.harmonics.constants import zeta3
+from ekore.harmonics.constants import zeta3
 
 
 def _flav_test(function):
diff --git a/tests/ekore/anomalous_dimensions/__init__.py b/tests/ekore/anomalous_dimensions/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/ekore/anomalous_dimensions/unpolarized/__init__.py b/tests/ekore/anomalous_dimensions/unpolarized/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py b/tests/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/eko/anomalous_dimensions/test_aem1.py b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_aem1.py
similarity index 57%
rename from tests/eko/anomalous_dimensions/test_aem1.py
rename to tests/ekore/anomalous_dimensions/unpolarized/space_like/test_aem1.py
index 70ee30cc5..0cf4d26e1 100644
--- a/tests/eko/anomalous_dimensions/test_aem1.py
+++ b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_aem1.py
@@ -1,22 +1,23 @@
 # Test LO splitting functions
 import numpy as np
 
-from eko import anomalous_dimensions as ad
+import ekore.harmonics as h
+import ekore.anomalous_dimensions.unpolarized.space_like as ad_us
 
 
 def test_number_conservation():
     # number
     N = complex(1.0, 0.0)
-    s1 = ad.harmonics.S1(N)
-    np.testing.assert_almost_equal(ad.aem1.gamma_ns(N, s1), 0)
+    s1 = h.S1(N)
+    np.testing.assert_almost_equal(ad_us.aem1.gamma_ns(N, s1), 0)
 
 
 def test_quark_momentum_conservation():
     # quark momentum
     N = complex(2.0, 0.0)
-    s1 = ad.harmonics.S1(N)
+    s1 = h.S1(N)
     np.testing.assert_almost_equal(
-        ad.aem1.gamma_ns(N, s1) + ad.aem1.gamma_phq(N),
+        ad_us.aem1.gamma_ns(N, s1) + ad_us.aem1.gamma_phq(N),
         0,
     )
 
@@ -26,5 +27,5 @@ def test_photon_momentum_conservation():
     N = complex(2.0, 0.0)
     for NF in range(2, 6 + 1):
         np.testing.assert_almost_equal(
-            ad.aem1.gamma_qph(N, NF) + ad.aem1.gamma_phph(NF), 0
+            ad_us.aem1.gamma_qph(N, NF) + ad_us.aem1.gamma_phph(NF), 0
         )
diff --git a/tests/eko/anomalous_dimensions/test_aem2.py b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_aem2.py
similarity index 93%
rename from tests/eko/anomalous_dimensions/test_aem2.py
rename to tests/ekore/anomalous_dimensions/unpolarized/space_like/test_aem2.py
index aadbbfc2f..0a1ddabde 100644
--- a/tests/eko/anomalous_dimensions/test_aem2.py
+++ b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_aem2.py
@@ -1,9 +1,9 @@
 # Test O(as1aem1) splitting functions
 import numpy as np
 
-from eko import anomalous_dimensions as ad
+import ekore.anomalous_dimensions.unpolarized.space_like as ad
 from eko import constants
-from eko import harmonics as h
+from ekore import harmonics as h
 
 
 def test_number_conservation():
diff --git a/tests/eko/anomalous_dimensions/test_as1.py b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as1.py
similarity index 90%
rename from tests/eko/anomalous_dimensions/test_as1.py
rename to tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as1.py
index 14fd3d004..23cd67d50 100644
--- a/tests/eko/anomalous_dimensions/test_as1.py
+++ b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as1.py
@@ -1,9 +1,8 @@
 # Test LO splitting functions
 import numpy as np
 
-import eko.anomalous_dimensions.aem1 as ad_aem1
-import eko.anomalous_dimensions.as1 as ad_as1
-from eko import harmonics
+import ekore.anomalous_dimensions.unpolarized.space_like.as1 as ad_as1
+from ekore import harmonics
 
 NF = 5
 
diff --git a/tests/eko/anomalous_dimensions/test_as1aem1.py b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as1aem1.py
similarity index 94%
rename from tests/eko/anomalous_dimensions/test_as1aem1.py
rename to tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as1aem1.py
index 988bb11f1..ce77d4f3a 100644
--- a/tests/eko/anomalous_dimensions/test_as1aem1.py
+++ b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as1aem1.py
@@ -2,9 +2,9 @@
 import numpy as np
 import pytest
 
-from eko import anomalous_dimensions as ad
+import ekore.anomalous_dimensions.unpolarized.space_like as ad
 from eko import constants
-from eko import harmonics as h
+from ekore import harmonics as h
 
 
 def test_number_conservation():
diff --git a/tests/eko/anomalous_dimensions/test_as2.py b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as2.py
similarity index 97%
rename from tests/eko/anomalous_dimensions/test_as2.py
rename to tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as2.py
index 01d7f2f28..a118e0a9a 100644
--- a/tests/eko/anomalous_dimensions/test_as2.py
+++ b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as2.py
@@ -1,8 +1,8 @@
 # Test NLO anomalous dims
 import numpy as np
 
-import eko.anomalous_dimensions.as2 as ad_as2
-import eko.harmonics as h
+import ekore.anomalous_dimensions.unpolarized.space_like.as2 as ad_as2
+import ekore.harmonics as h
 from eko import constants as const
 
 NF = 5
diff --git a/tests/eko/anomalous_dimensions/test_as3.py b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as3.py
similarity index 89%
rename from tests/eko/anomalous_dimensions/test_as3.py
rename to tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as3.py
index e0830c4a2..7db2b7c64 100644
--- a/tests/eko/anomalous_dimensions/test_as3.py
+++ b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as3.py
@@ -1,8 +1,8 @@
 # Test NNLO anomalous dimensions
 import numpy as np
 
-import eko.anomalous_dimensions.as3 as ad_as3
-from eko import harmonics as h
+import ekore.anomalous_dimensions.unpolarized.space_like.as3 as ad_as3
+from ekore import harmonics as h
 
 NF = 5
 
diff --git a/tests/eko/anomalous_dimensions/test_as4.py b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as4.py
similarity index 98%
rename from tests/eko/anomalous_dimensions/test_as4.py
rename to tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as4.py
index 95c02f180..38935ec00 100644
--- a/tests/eko/anomalous_dimensions/test_as4.py
+++ b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as4.py
@@ -1,7 +1,7 @@
 # Test N3LO anomalous dimensions
 import numpy as np
 
-from eko.anomalous_dimensions.as4 import (
+from ekore.anomalous_dimensions.unpolarized.space_like.as4 import (
     gamma_singlet,
     ggg,
     ggq,
@@ -12,7 +12,7 @@
     gqg,
 )
 from eko.constants import CA, CF
-from eko.harmonics import compute_cache
+from ekore.harmonics import compute_cache
 
 NF = 5
 
diff --git a/tests/eko/anomalous_dimensions/test_init.py b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py
similarity index 74%
rename from tests/eko/anomalous_dimensions/test_init.py
rename to tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py
index a080cdd1b..f36efec04 100644
--- a/tests/eko/anomalous_dimensions/test_init.py
+++ b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py
@@ -4,10 +4,11 @@
 import numpy as np
 from numpy.testing import assert_allclose, assert_almost_equal, assert_raises
 
-from eko import anomalous_dimensions as ad
+from ekore import anomalous_dimensions as ad
 from eko import basis_rotation as br
-from eko.anomalous_dimensions import as1 as ad_as1
-from eko.anomalous_dimensions import harmonics
+import ekore.anomalous_dimensions.unpolarized.space_like as ad_us
+from ekore.anomalous_dimensions.unpolarized.space_like import as1 as ad_as1
+from ekore import harmonics
 
 NF = 5
 
@@ -41,7 +42,7 @@ def test_eigensystem_gamma_singlet_projectors_EV():
             # NNLO and N3LO too big numbers,
             # ignore Runtime Warnings
             warnings.simplefilter("ignore", RuntimeWarning)
-            for gamma_S in ad.gamma_singlet(o, N, nf):
+            for gamma_S in ad_us.gamma_singlet(o, N, nf):
                 _exp, l_p, l_m, e_p, e_m = ad.exp_singlet(gamma_S)
                 # projectors behave as P_a . P_b = delta_ab P_a
                 assert_allclose(np.dot(e_p, e_p), e_p)
@@ -56,40 +57,40 @@ def test_gamma_ns():
     nf = 3
     # LO
     assert_almost_equal(
-        ad.gamma_ns((3, 0), br.non_singlet_pids_map["ns+"], 1, nf)[0], 0.0
+        ad_us.gamma_ns((3, 0), br.non_singlet_pids_map["ns+"], 1, nf)[0], 0.0
     )
     # NLO
     assert_allclose(
-        ad.gamma_ns((2, 0), br.non_singlet_pids_map["ns-"], 1, nf),
+        ad_us.gamma_ns((2, 0), br.non_singlet_pids_map["ns-"], 1, nf),
         np.zeros(2),
         atol=2e-6,
     )
     # NNLO
     assert_allclose(
-        ad.gamma_ns((3, 0), br.non_singlet_pids_map["ns-"], 1, nf),
+        ad_us.gamma_ns((3, 0), br.non_singlet_pids_map["ns-"], 1, nf),
         np.zeros(3),
         atol=2e-4,
     )
     assert_allclose(
-        ad.gamma_ns((3, 0), br.non_singlet_pids_map["nsV"], 1, nf),
+        ad_us.gamma_ns((3, 0), br.non_singlet_pids_map["nsV"], 1, nf),
         np.zeros(3),
         atol=8e-4,
     )
     # N3LO
     assert_allclose(
-        ad.gamma_ns((4, 0), br.non_singlet_pids_map["ns-"], 1, nf),
+        ad_us.gamma_ns((4, 0), br.non_singlet_pids_map["ns-"], 1, nf),
         np.zeros(4),
         atol=2e-4,
     )
     # N3LO valence has a spurious pole, need to add a small shift
     assert_allclose(
-        ad.gamma_ns((4, 0), br.non_singlet_pids_map["nsV"], 1 + 1e-6, nf),
+        ad_us.gamma_ns((4, 0), br.non_singlet_pids_map["nsV"], 1 + 1e-6, nf),
         np.zeros(4),
         atol=5e-4,
     )
     assert_raises(
         AssertionError,
         assert_allclose,
-        ad.gamma_ns((4, 0), br.non_singlet_pids_map["ns+"], 1, nf),
+        ad_us.gamma_ns((4, 0), br.non_singlet_pids_map["ns+"], 1, nf),
         np.zeros(4),
     )
diff --git a/tests/ekore/harmonics/__init__.py b/tests/ekore/harmonics/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/eko/harmonics/test_f_functions.py b/tests/ekore/harmonics/test_f_functions.py
similarity index 98%
rename from tests/eko/harmonics/test_f_functions.py
rename to tests/ekore/harmonics/test_f_functions.py
index 943dd7131..89136a9a6 100644
--- a/tests/eko/harmonics/test_f_functions.py
+++ b/tests/ekore/harmonics/test_f_functions.py
@@ -2,8 +2,8 @@
 
 import numpy as np
 
-from eko import harmonics
-from eko.harmonics import w5
+from ekore import harmonics
+from ekore.harmonics import w5
 
 zeta2 = harmonics.constants.zeta2
 zeta3 = harmonics.constants.zeta3
diff --git a/tests/eko/harmonics/test_g_functions.py b/tests/ekore/harmonics/test_g_functions.py
similarity index 98%
rename from tests/eko/harmonics/test_g_functions.py
rename to tests/ekore/harmonics/test_g_functions.py
index 4dc1e758f..b6a895e6e 100644
--- a/tests/eko/harmonics/test_g_functions.py
+++ b/tests/ekore/harmonics/test_g_functions.py
@@ -2,7 +2,7 @@
 
 import numpy as np
 
-from eko.anomalous_dimensions import harmonics as h
+from ekore import harmonics as h
 
 zeta3 = h.constants.zeta3
 log2 = h.constants.log2
diff --git a/tests/eko/harmonics/test_init.py b/tests/ekore/harmonics/test_init.py
similarity index 99%
rename from tests/eko/harmonics/test_init.py
rename to tests/ekore/harmonics/test_init.py
index fbaa4cea4..dc4e39637 100644
--- a/tests/eko/harmonics/test_init.py
+++ b/tests/ekore/harmonics/test_init.py
@@ -1,6 +1,6 @@
 import numpy as np
 
-from eko import harmonics as h
+from ekore import harmonics as h
 
 
 def test_spm1():
diff --git a/tests/eko/harmonics/test_log_functions.py b/tests/ekore/harmonics/test_log_functions.py
similarity index 97%
rename from tests/eko/harmonics/test_log_functions.py
rename to tests/ekore/harmonics/test_log_functions.py
index e53f13c11..2f3ad4e90 100644
--- a/tests/eko/harmonics/test_log_functions.py
+++ b/tests/ekore/harmonics/test_log_functions.py
@@ -1,7 +1,7 @@
 import numpy as np
 from scipy.integrate import quad
 
-from eko import harmonics as h
+from ekore import harmonics as h
 
 
 def test_lm1pm1():
diff --git a/tests/eko/harmonics/test_polygamma.py b/tests/ekore/harmonics/test_polygamma.py
similarity index 99%
rename from tests/eko/harmonics/test_polygamma.py
rename to tests/ekore/harmonics/test_polygamma.py
index a790d3828..7ea4389dd 100644
--- a/tests/eko/harmonics/test_polygamma.py
+++ b/tests/ekore/harmonics/test_polygamma.py
@@ -1,7 +1,7 @@
 import numpy as np
 import pytest
 
-from eko import harmonics
+from ekore import harmonics
 
 # until https://github.com/numba/numba/pull/5660 is confirmed
 # we need to deactivate numba prior running
diff --git a/tests/ekore/operator_matrix_elements/__init__.py b/tests/ekore/operator_matrix_elements/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/ekore/operator_matrix_elements/unpolarized/__init__.py b/tests/ekore/operator_matrix_elements/unpolarized/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/ekore/operator_matrix_elements/unpolarized/space_like/__init__.py b/tests/ekore/operator_matrix_elements/unpolarized/space_like/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/eko/matching_conditions/test_n3lo.py b/tests/ekore/operator_matrix_elements/unpolarized/space_like/test_n3lo.py
similarity index 97%
rename from tests/eko/matching_conditions/test_n3lo.py
rename to tests/ekore/operator_matrix_elements/unpolarized/space_like/test_n3lo.py
index 377e1b606..eed992f3d 100644
--- a/tests/eko/matching_conditions/test_n3lo.py
+++ b/tests/ekore/operator_matrix_elements/unpolarized/space_like/test_n3lo.py
@@ -1,9 +1,9 @@
 # Test N3LO OME
 import numpy as np
 
-from eko.harmonics import compute_cache
-from eko.matching_conditions import as3
-from eko.matching_conditions.as3 import A_ns, A_qqNS, A_singlet
+from ekore.harmonics import compute_cache
+from ekore.operator_matrix_elements.unpolarized.space_like import as3
+from ekore.operator_matrix_elements.unpolarized.space_like.as3 import A_ns, A_qqNS, A_singlet
 
 
 def test_A_3():
diff --git a/tests/eko/matching_conditions/test_nlo.py b/tests/ekore/operator_matrix_elements/unpolarized/space_like/test_nlo.py
similarity index 91%
rename from tests/eko/matching_conditions/test_nlo.py
rename to tests/ekore/operator_matrix_elements/unpolarized/space_like/test_nlo.py
index b7f03cc82..953d797c7 100644
--- a/tests/eko/matching_conditions/test_nlo.py
+++ b/tests/ekore/operator_matrix_elements/unpolarized/space_like/test_nlo.py
@@ -1,8 +1,8 @@
 # Test NLO OME
 import numpy as np
 
-from eko.harmonics import compute_cache
-from eko.matching_conditions.as1 import A_ns, A_singlet
+from ekore.harmonics import compute_cache
+from ekore.operator_matrix_elements.unpolarized.space_like.as1 import A_ns, A_singlet
 
 
 def test_A_1_intrinsic():
diff --git a/tests/eko/matching_conditions/test_nnlo.py b/tests/ekore/operator_matrix_elements/unpolarized/space_like/test_nnlo.py
similarity index 96%
rename from tests/eko/matching_conditions/test_nnlo.py
rename to tests/ekore/operator_matrix_elements/unpolarized/space_like/test_nnlo.py
index aea957e95..648568d1c 100644
--- a/tests/eko/matching_conditions/test_nnlo.py
+++ b/tests/ekore/operator_matrix_elements/unpolarized/space_like/test_nnlo.py
@@ -2,8 +2,8 @@
 
 import numpy as np
 
-from eko.harmonics import compute_cache, constants
-from eko.matching_conditions.as2 import A_ns, A_qq_ns, A_singlet
+from ekore.harmonics import compute_cache, constants
+from ekore.operator_matrix_elements.unpolarized.space_like.as2 import A_ns, A_qq_ns, A_singlet
 
 
 def test_A_2():