diff --git a/benchmarks/eko/benchmark_evol_to_unity.py b/benchmarks/eko/benchmark_evol_to_unity.py
index 6baf6bfa8..e780fce6d 100644
--- a/benchmarks/eko/benchmark_evol_to_unity.py
+++ b/benchmarks/eko/benchmark_evol_to_unity.py
@@ -55,8 +55,9 @@ def test_operator_grid(
         q20 = 30
         q21 = 50
         nf = 4
-        o = Operator(g.config, g.managers, nf, q20, q21)
-        o_back = Operator(g.config, g.managers, nf, q21, q20)
+        p = False
+        o = Operator(g.config, g.managers, nf, p, q20, q21)
+        o_back = Operator(g.config, g.managers, nf, p, q21, q20)
         o.compute()
         o_back.compute()
 
diff --git a/benchmarks/performance/harmonics.py b/benchmarks/performance/harmonics.py
index f803bcda0..e8df08b87 100644
--- a/benchmarks/performance/harmonics.py
+++ b/benchmarks/performance/harmonics.py
@@ -17,12 +17,12 @@ def setup(self):
 
     def time_as1_sing(self):
         for n in self.ns:
-            ad.gamma_singlet(0, n, NF)
+            ad.gamma_singlet(0, n, NF, p=False)
 
     def time_as2_sing(self):
         for n in self.ns:
-            ad.gamma_singlet(1, n, NF)
+            ad.gamma_singlet(1, n, NF, p=False)
 
     def time_as3_sing(self):
         for n in self.ns:
-            ad.gamma_singlet(2, n, NF)
+            ad.gamma_singlet(2, n, NF, p=False)
diff --git a/src/eko/anomalous_dimensions/__init__.py b/src/eko/anomalous_dimensions/__init__.py
index c3a3f8c78..6150ea7d5 100644
--- a/src/eko/anomalous_dimensions/__init__.py
+++ b/src/eko/anomalous_dimensions/__init__.py
@@ -21,7 +21,7 @@
 
 from .. import basis_rotation as br
 from .. import harmonics
-from . import aem1, aem2, as1, as1aem1, as2, as3, as4
+from . import aem1, aem2, as1, as1aem1, as2, as3, as4, asp1
 
 
 @nb.njit(cache=True)
@@ -73,7 +73,7 @@ def exp_singlet(gamma_S):
 
 
 @nb.njit(cache=True)
-def gamma_ns(order, mode, n, nf):
+def gamma_ns(order, mode, n, nf, p):
     r"""Computes the tower of the non-singlet anomalous dimensions
 
     Parameters
@@ -86,6 +86,8 @@ def gamma_ns(order, mode, n, nf):
         Mellin variable
     nf : int
         Number of active flavors
+    p : Boolean
+        Polarised (True) or un-Polarised (False)
 
     Returns
     -------
@@ -120,9 +122,12 @@ def gamma_ns(order, mode, n, nf):
         sx = harmonics.sx(n, max_weight=order[0] + 1)
     # now combine
     gamma_ns = np.zeros(order[0], np.complex_)
-    gamma_ns[0] = as1.gamma_ns(n, sx[0])
-    # NLO and beyond
-    if order[0] >= 2:
+    if p == False:
+        gamma_ns[0] = as1.gamma_ns(n, sx[0])
+    elif p == True:
+        gamma_ns[0] = asp1.gamma_pns(n, sx[0])
+    # NLO and beyondd
+    if order[0] >= 2 and p == False:
         if mode == 10101:
             gamma_ns_1 = as2.gamma_nsp(n, nf, sx)
         # To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here
@@ -131,8 +136,12 @@ def gamma_ns(order, mode, n, nf):
         else:
             raise NotImplementedError("Non-singlet sector is not implemented")
         gamma_ns[1] = gamma_ns_1
+
+    elif order[0] >= 2 and p == True:
+        raise Exception("Polarised case is not known at this order")
+
     # NNLO and beyond
-    if order[0] >= 3:
+    if order[0] >= 3 and p == False:
         if mode == 10101:
             gamma_ns_2 = as3.gamma_nsp(n, nf, sx)
         elif mode == 10201:
@@ -140,8 +149,12 @@ def gamma_ns(order, mode, n, nf):
         elif mode == 10200:
             gamma_ns_2 = as3.gamma_nsv(n, nf, sx)
         gamma_ns[2] = gamma_ns_2
+
+    elif order[0] >= 3 and p == True:
+        raise Exception("Polarised case is not known at this order")
+
     # N3LO
-    if order[0] >= 4:
+    if order[0] >= 4 and p == False:
         if mode == 10101:
             gamma_ns_3 = as4.gamma_nsp(n, nf, full_sx_cache)
         elif mode == 10201:
@@ -149,11 +162,15 @@ def gamma_ns(order, mode, n, nf):
         elif mode == 10200:
             gamma_ns_3 = as4.gamma_nsv(n, nf, full_sx_cache)
         gamma_ns[3] = gamma_ns_3
+
+    elif order[0] >= 4 and p == True:
+        raise Exception("Polarised case is not known at this order")
+
     return gamma_ns
 
 
 @nb.njit(cache=True)
-def gamma_singlet(order, n, nf):
+def gamma_singlet(order, n, nf, p):
     r"""Computes the tower of the singlet anomalous dimensions matrices
 
     Parameters
@@ -164,7 +181,8 @@ def gamma_singlet(order, n, nf):
         Mellin variable
     nf : int
         Number of active flavors
-
+    p : Boolean
+        Polarised (True) or un-Polarised (False)
     Returns
     -------
     numpy.ndarray
@@ -196,11 +214,25 @@ def gamma_singlet(order, n, nf):
         sx = harmonics.sx(n, max_weight=order[0])
 
     gamma_s = np.zeros((order[0], 2, 2), np.complex_)
-    gamma_s[0] = as1.gamma_singlet(n, sx[0], nf)
-    if order[0] >= 2:
+
+    if p == False:
+        gamma_s[0] = as1.gamma_singlet(n, sx[0], nf)
+    elif p == True:
+        gamma_s[0] = asp1.gamma_psinglet(n, sx[0], nf)
+
+    if order[0] >= 2 and p == False:
         gamma_s[1] = as2.gamma_singlet(n, nf, sx)
-    if order[0] >= 3:
+    elif order[0] >= 2 and p == True:
+        raise Exception("Polarised case is not known at this order")
+
+    if order[0] >= 3 and p == False:
         gamma_s[2] = as3.gamma_singlet(n, nf, sx)
-    if order[0] >= 4:
+    elif order[0] >= 3 and p == True:
+        raise Exception("Polarised case is not known at this order")
+
+    if order[0] >= 4 and p == False:
         gamma_s[3] = as4.gamma_singlet(n, nf, full_sx_cache)
+    elif order[0] >= 4 and p == True:
+        raise Exception("Polarised case is not known at this order")
+
     return gamma_s
diff --git a/src/eko/anomalous_dimensions/aem1.py b/src/eko/anomalous_dimensions/aem1.py
index f272439cb..36fbb5c6b 100644
--- a/src/eko/anomalous_dimensions/aem1.py
+++ b/src/eko/anomalous_dimensions/aem1.py
@@ -8,6 +8,7 @@
 from . import as1
 
 
+
 @nb.njit(cache=True)
 def gamma_phq(N):
     """
diff --git a/src/eko/anomalous_dimensions/asp1.py b/src/eko/anomalous_dimensions/asp1.py
new file mode 100644
index 000000000..030fda480
--- /dev/null
+++ b/src/eko/anomalous_dimensions/asp1.py
@@ -0,0 +1,127 @@
+# -*- coding: utf-8 -*-
+"""This file contains the leading-order polarised  Altarelli-Parisi splitting kernels."""
+
+import numba as nb
+import numpy as np
+
+from .. import constants
+from .as1 import gamma_ns as gamma_pns
+
+# @nb.njit(cache=True)
+# def gamma_pns(N, s1):
+# Computes the leading-order non-singlet anomalous dimension for the polarised case.
+# This is going to be the same expression as the one for the unpolarised case.
+
+# gamma = -(3.0 - 4.0 * s1 + 2.0 / N / (N + 1.0))
+# result = constants.CF * gamma
+# return result
+
+
+@nb.njit(cache=True)
+def gamma_pqg(N, nf):
+    """
+    Computes the leading-order polarised quark-gluon anomalous dimension
+    #requires citation
+
+    Parameters
+    ----------
+      N : complex
+        Mellin moment
+      nf : int
+        Number of active flavors
+
+    Returns
+    -------
+      gamma_qg : complex
+        Leading-order polarised quark-gluon anomalous dimension :math:`\\gamma_{qg}^{(0)}(N)`
+    """
+    gamma = -(N - 1) / N / (N + 1)
+    result = 2.0 * constants.TR * 2.0 * nf * gamma
+    return result
+
+
+@nb.njit(cache=True)
+def gamma_pgq(N):
+    """
+    Computes the leading-order polarised gluon-quark anomalous dimension
+
+
+    Parameters
+    ----------
+      N : complex
+        Mellin moment
+
+    Returns
+    -------
+      gamma_gq : complex
+        Leading-order gluon-quark anomalous dimension :math:`\\gamma_{gq}^{(0)}(N)`
+    """
+    gamma = -(N + 2) / N / (N + 1)
+    result = 2.0 * constants.CF * gamma
+    return result
+
+
+@nb.njit(cache=True)
+def gamma_pgg(N, s1, nf):
+    """
+    Computes the leading-order polarised gluon-gluon anomalous dimension
+
+
+    Parameters
+    ----------
+      N : complex
+        Mellin moment
+      s1 : complex
+        harmonic sum :math:`S_{1}`
+      nf : int
+        Number of active flavors
+
+    Returns
+    -------
+      gamma_gg : complex
+        Leading-order gluon-gluon anomalous dimension :math:`\\gamma_{gg}^{(0)}(N)`
+    """
+    gamma = s1 - 2 / N / (N + 1)
+    result = constants.CA * (4.0 * gamma - 11.0 / 3.0) + 4.0 / 3.0 * constants.TR * nf
+    return result
+
+
+
+@nb.njit(cache=True)
+def gamma_psinglet(N, s1, nf):
+    r"""
+      Computes the leading-order polarised singlet anomalous dimension matrix
+
+      .. math::
+          \gamma_S^{(0)} = \left(\begin{array}{cc}
+            \gamma_{qq}^{(0)} & \gamma_{qg}^{(0)}\\
+            \gamma_{gq}^{(0)} & \gamma_{gg}^{(0)}
+          \end{array}\right)
+
+      Parameters
+      ----------
+        N : complex
+          Mellin moment
+        s1 : complex
+          harmonic sum :math:`S_{1}`
+        nf : int
+          Number of active flavors
+
+      Returns
+      -------
+        gamma_S_0 : numpy.ndarray
+          Leading-order singlet anomalous dimension matrix :math:`\gamma_{S}^{(0)}(N)`
+
+      See Also
+      --------
+        gamma_ns : :math:`\gamma_{qq}^{(0)}`
+        gamma_qg : :math:`\gamma_{qg}^{(0)}`
+        gamma_gq : :math:`\gamma_{gq}^{(0)}`
+        gamma_gg : :math:`\gamma_{gg}^{(0)}`
+    """
+    gamma_pqq = gamma_pns(N, s1)
+    gamma_pS_0 = np.array(
+        [[gamma_pqq, gamma_pqg(N, nf)], [gamma_pgq(N), gamma_pgg(N, s1, nf)]],
+        np.complex_,
+    )
+    return gamma_pS_0
diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index 323b117b4..644c96564 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -11,8 +11,6 @@
 
 import numba as nb
 import numpy as np
-from scipy import integrate
-
 from .. import anomalous_dimensions as ad
 from .. import basis_rotation as br
 from .. import interpolation, mellin
@@ -20,6 +18,7 @@
 from ..kernels import non_singlet as ns
 from ..kernels import singlet as s
 from ..member import OpMember
+from scipy import integrate
 
 logger = logging.getLogger(__name__)
 
@@ -84,7 +83,7 @@ def path(self):
 
     @property
     def n(self):
-        """Returs the Mellin moment :math:`N`."""
+        """Return the Mellin moment :math:`N`."""
         return self.path.n
 
     def integrand(
@@ -125,6 +124,7 @@ def quad_ker(
     as0,
     as_raw,
     nf,
+    p,
     L,
     ev_op_iterations,
     ev_op_max_order,
@@ -159,6 +159,8 @@ def quad_ker(
         coupling value at the process scale
     nf : int
         number of active flavors
+    p : Boolean
+        Polarised (True) or un-Polarised (False)
     L : float
         logarithm of the squared ratio of factorization and renormalization scale
     ev_op_iterations : int
@@ -182,7 +184,7 @@ def quad_ker(
 
     # compute the actual evolution kernel
     if ker_base.is_singlet:
-        gamma_singlet = ad.gamma_singlet(order, ker_base.n, nf)
+        gamma_singlet = ad.gamma_singlet(order, ker_base.n, nf, p)
         # scale var exponentiated is directly applied on gamma
         if sv_mode == sv.Modes.exponentiated:
             gamma_singlet = sv.exponentiated.gamma_variation(
@@ -198,14 +200,14 @@ def quad_ker(
             ev_op_iterations,
             ev_op_max_order,
         )
-        # scale var expanded is applied on the kernel
+
         if sv_mode == sv.Modes.expanded and not is_threshold:
             ker = np.ascontiguousarray(
                 sv.expanded.singlet_variation(gamma_singlet, as_raw, order, nf, L)
             ) @ np.ascontiguousarray(ker)
         ker = select_singlet_element(ker, mode0, mode1)
     else:
-        gamma_ns = ad.gamma_ns(order, mode0, ker_base.n, nf)
+        gamma_ns = ad.gamma_ns(order, mode0, ker_base.n, nf, p)
         if sv_mode == sv.Modes.exponentiated:
             gamma_ns = sv.exponentiated.gamma_variation(gamma_ns, order, nf, L)
         ker = ns.dispatcher(
@@ -254,7 +256,14 @@ class Operator(sv.ModeMixin):
     full_labels = br.full_labels
 
     def __init__(
-        self, config, managers, nf, q2_from, q2_to, mellin_cut=5e-2, is_threshold=False
+        self,
+        config,
+        managers,
+        nf,
+        q2_from,
+        q2_to,
+        mellin_cut=5e-2,
+        is_threshold=False,
     ):
         self.config = config
         self.managers = managers
@@ -377,7 +386,8 @@ def quad_ker(self, label, logx, areas):
             as0=self.a_s[0],
             as_raw=self.a_s[2],
             nf=self.nf,
-            L=np.log(self.xif2),
+            p=self.config["p"],
+            L=np.log(self.fact_to_ren),
             ev_op_iterations=self.config["ev_op_iterations"],
             ev_op_max_order=tuple(self.config["ev_op_max_order"]),
             sv_mode=self.sv_mode,
diff --git a/src/eko/evolution_operator/grid.py b/src/eko/evolution_operator/grid.py
index 1b55a2072..ea2b7fc08 100644
--- a/src/eko/evolution_operator/grid.py
+++ b/src/eko/evolution_operator/grid.py
@@ -7,7 +7,6 @@
 
 import logging
 import numbers
-
 import numpy as np
 import numpy.typing as npt
 
@@ -82,6 +81,7 @@ def __init__(
 
         for i, q in enumerate("cbt"):
             config[f"m{q}"] = masses[i]
+        p = config["p"]
         method = config["method"] = configs.evolution_method.value
         config["backward_inversion"] = configs.inversion_method
         config["ev_op_max_order"] = configs.ev_op_max_order
@@ -223,7 +223,11 @@ def generate(self, q2):
         thr_ops = self.get_threshold_operators(path)
         # we start composing with the highest operator ...
         operator = Operator(
-            self.config, self.managers, path[-1].nf, path[-1].q2_from, path[-1].q2_to
+            self.config,
+            self.managers,
+            path[-1].nf,
+            path[-1].q2_from,
+            path[-1].q2_to,
         )
         operator.compute()
         intrinsic_range = self.config["intrinsic_range"]
diff --git a/tests/eko/anomalous_dimensions/test_init.py b/tests/eko/anomalous_dimensions/test_init.py
index a080cdd1b..e8c768bfc 100644
--- a/tests/eko/anomalous_dimensions/test_init.py
+++ b/tests/eko/anomalous_dimensions/test_init.py
@@ -36,12 +36,13 @@ def test_eigensystem_gamma_singlet_0_values():
 
 def test_eigensystem_gamma_singlet_projectors_EV():
     nf = 3
+    p = False
     for N in [3, 4]:  # N=2 seems close to 0, so test fails
         for o in [(2, 0), (3, 0), (4, 0)]:
             # 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.gamma_singlet(o, N, nf, p):
                 _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)
@@ -54,42 +55,44 @@ def test_eigensystem_gamma_singlet_projectors_EV():
 
 def test_gamma_ns():
     nf = 3
+    p = False
     # LO
     assert_almost_equal(
-        ad.gamma_ns((3, 0), br.non_singlet_pids_map["ns+"], 1, nf)[0], 0.0
+        ad.gamma_ns((3, 0), br.non_singlet_pids_map["ns+"], 1, nf, p)[0],
+        0.0,
     )
     # NLO
     assert_allclose(
-        ad.gamma_ns((2, 0), br.non_singlet_pids_map["ns-"], 1, nf),
+        ad.gamma_ns((2, 0), br.non_singlet_pids_map["ns-"], 1, nf, p),
         np.zeros(2),
         atol=2e-6,
     )
     # NNLO
     assert_allclose(
-        ad.gamma_ns((3, 0), br.non_singlet_pids_map["ns-"], 1, nf),
+        ad.gamma_ns((3, 0), br.non_singlet_pids_map["ns-"], 1, nf, p),
         np.zeros(3),
         atol=2e-4,
     )
     assert_allclose(
-        ad.gamma_ns((3, 0), br.non_singlet_pids_map["nsV"], 1, nf),
+        ad.gamma_ns((3, 0), br.non_singlet_pids_map["nsV"], 1, nf, p),
         np.zeros(3),
         atol=8e-4,
     )
     # N3LO
     assert_allclose(
-        ad.gamma_ns((4, 0), br.non_singlet_pids_map["ns-"], 1, nf),
+        ad.gamma_ns((4, 0), br.non_singlet_pids_map["ns-"], 1, nf, p),
         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.gamma_ns((4, 0), br.non_singlet_pids_map["nsV"], 1 + 1e-6, nf, p),
         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.gamma_ns((4, 0), br.non_singlet_pids_map["ns+"], 1, nf, p),
         np.zeros(4),
     )
diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py
index dc8ddbd8e..27f219d72 100644
--- a/tests/eko/evolution_operator/test_init.py
+++ b/tests/eko/evolution_operator/test_init.py
@@ -39,6 +39,7 @@ def test_quad_ker(monkeypatch):
             as0=2,
             as_raw=1,
             nf=3,
+            p=False,
             L=0,
             ev_op_iterations=0,
             ev_op_max_order=(0, 0),
@@ -59,6 +60,7 @@ def test_quad_ker(monkeypatch):
             as0=2,
             as_raw=1,
             nf=3,
+            p=False,
             L=0,
             ev_op_iterations=0,
             ev_op_max_order=(0, 0),
@@ -79,6 +81,7 @@ def test_quad_ker(monkeypatch):
             as0=2,
             as_raw=1,
             nf=3,
+            p=False,
             L=0,
             ev_op_iterations=0,
             ev_op_max_order=(0, 0),
@@ -101,6 +104,7 @@ def test_quad_ker(monkeypatch):
                 as0=2,
                 as_raw=1,
                 nf=3,
+                p=False,
                 L=0,
                 ev_op_iterations=0,
                 ev_op_max_order=(1, 0),
@@ -123,6 +127,7 @@ def test_quad_ker(monkeypatch):
         as0=2,
         as_raw=1,
         nf=3,
+        p=False,
         L=0,
         ev_op_iterations=0,
         ev_op_max_order=(0, 0),
@@ -286,13 +291,13 @@ def test_compute(self, monkeypatch, theory_ffns, operator_card, tmp_path):
 
 def test_pegasus_path():
     def quad_ker_pegasus(
-        u, order, mode0, method, logx, areas, a1, a0, nf, ev_op_iterations
+        u, order, mode0, method, logx, areas, a1, a0, nf, p, ev_op_iterations
     ):
         # compute the mellin inversion as done in pegasus
         phi = 3 / 4 * np.pi
         c = 1.9
         n = complex(c + u * np.exp(1j * phi))
-        gamma_ns = ad.gamma_ns(order, mode0, n, nf)
+        gamma_ns = ad.gamma_ns(order, mode0, n, nf, p)
         ker = ns.dispatcher(
             order,
             method,
@@ -317,6 +322,7 @@ def quad_ker_pegasus(
     as_raw = a1 = 1
     a0 = 2
     nf = 3
+    p = False
     L = 0
     ev_op_iterations = 10
     for logx in logxs:
@@ -337,6 +343,7 @@ def quad_ker_pegasus(
                     a0,
                     as_raw,
                     nf,
+                    p,
                     L,
                     ev_op_iterations,
                     10,
@@ -362,6 +369,7 @@ def quad_ker_pegasus(
                     a1,
                     a0,
                     nf,
+                    p,
                     ev_op_iterations,
                 ),
                 epsabs=1e-12,
diff --git a/tests/eko/kernels/test_ns.py b/tests/eko/kernels/test_ns.py
index 17d4964b8..d3ad36908 100644
--- a/tests/eko/kernels/test_ns.py
+++ b/tests/eko/kernels/test_ns.py
@@ -206,7 +206,7 @@ def test_error():
     with pytest.raises(NotImplementedError):
         ns.dispatcher((5, 0), "iterate-exact", np.random.rand(3) + 0j, 0.2, 0.1, 3, 10)
     with pytest.raises(NotImplementedError):
-        ad.gamma_ns((2, 0), 10202, 1, 3)
+        ad.gamma_ns((2, 0), 10202, 1, 3, p=False)
 
 
 def test_gamma_usage():
diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py
index bce253f9c..edda0db08 100644
--- a/tests/eko/scale_variations/test_expanded.py
+++ b/tests/eko/scale_variations/test_expanded.py
@@ -79,7 +79,7 @@ def scheme_diff(g, k, pto, is_singlet):
     for L in [np.log(0.5), np.log(2)]:
         for order in [(2, 0), (3, 0)]:
             # Non singlet kernels
-            gns = gamma_ns(order, br.non_singlet_pids_map["ns+"], n, nf)
+            gns = gamma_ns(order, br.non_singlet_pids_map["ns+"], n, nf, p=False)
             ker = non_singlet.dispatcher(
                 order, method, gns, a1, a0, nf, ev_op_iterations=1
             )
@@ -92,7 +92,7 @@ def scheme_diff(g, k, pto, is_singlet):
             np.testing.assert_allclose((ker_a - ker_b) / ker, ns_diff, atol=1e-3)
 
             # Singlet kernels
-            gs = gamma_singlet(order, n, nf)
+            gs = gamma_singlet(order, n, nf, p=False)
             ker = singlet.dispatcher(
                 order, method, gs, a1, a0, nf, ev_op_iterations=1, ev_op_max_order=1
             )
diff --git a/tests/eko/test_ad_asp1.py b/tests/eko/test_ad_asp1.py
new file mode 100755
index 000000000..09581a201
--- /dev/null
+++ b/tests/eko/test_ad_asp1.py
@@ -0,0 +1,25 @@
+# Test LO Polarised splitting functions
+import numpy as np
+import eko.anomalous_dimensions.asp1 as ad_asp1
+from eko import harmonics
+
+NF = 5
+
+
+def test_quark_momentum_conservation():
+    # quark momentum
+    N = complex(2.0, 0.0)
+    s1 = harmonics.S1(N)
+    np.testing.assert_almost_equal(
+        ad_asp1.gamma_pns(N, s1) + ad_asp1.gamma_pgq(N),
+        0,
+    )
+
+
+def test_gluon_momentum_conservation():
+    # gluon momentum
+    N = complex(2.0, 0.0)
+    s1 = harmonics.S1(N)
+    np.testing.assert_almost_equal(
+        ad_asp1.gamma_pqg(N, NF) + ad_asp1.gamma_pgg(N, s1, NF), 0
+    )