diff --git a/benchmarks/eko/benchmark_ad.py b/benchmarks/eko/benchmark_ad.py
index cb5c07e66..f8ba47bd2 100644
--- a/benchmarks/eko/benchmark_ad.py
+++ b/benchmarks/eko/benchmark_ad.py
@@ -138,8 +138,9 @@ def check_gamma_1_pegasus(N, NF):
     P1NSM = CF * ((CF - CA / 2.0) * PNMA + CA * PNSB + TR * NF * PNSC)
 
     sx = h.sx(N, 2)
-    np.testing.assert_allclose(ad_as2.gamma_nsp(N, NF, sx), -P1NSP)
-    np.testing.assert_allclose(ad_as2.gamma_nsm(N, NF, sx), -P1NSM)
+    sx_new = h.compute_additional_sx_cache(N, False, False, sx)
+    np.testing.assert_allclose(ad_as2.gamma_nsp(N, NF, sx, sx_new), -P1NSP)
+    np.testing.assert_allclose(ad_as2.gamma_nsm(N, NF, sx, sx_new), -P1NSM)
 
     NS = N * N
     NT = NS * N
@@ -256,7 +257,8 @@ def check_gamma_1_pegasus(N, NF):
     P1Sgq = (CF * CF * PGQA + CF * CA * PGQB + TR * NF * CF * PGQC) * 4.0
     P1Sgg = (CA * CA * PGGA + TR * NF * (CA * PGGB + CF * PGGC)) * 4.0
 
-    gS1 = ad_as2.gamma_singlet(N, NF, sx)
+    sx_new = h.compute_additional_sx_cache(N, True, False, sx)
+    gS1 = ad_as2.gamma_singlet(N, NF, sx, sx_new)
     np.testing.assert_allclose(gS1[0, 0], -P1Sqq)
     np.testing.assert_allclose(gS1[0, 1], -P1Sqg)
     np.testing.assert_allclose(gS1[1, 0], -P1Sgq)
diff --git a/src/eko/anomalous_dimensions/__init__.py b/src/eko/anomalous_dimensions/__init__.py
index bd5dae69c..bf0ee108d 100644
--- a/src/eko/anomalous_dimensions/__init__.py
+++ b/src/eko/anomalous_dimensions/__init__.py
@@ -135,7 +135,9 @@ def gamma_ns(order, mode, n, nf):
     """
     # cache the s-es
     if order[0] >= 4:
-        full_sx_cache = harmonics.compute_cache(n, 5, is_singlet=False)
+        full_sx_cache = harmonics.compute_cache(
+            n, 5, is_singlet=False, additional=True, qed=False
+        )
         sx = np.array(
             [
                 full_sx_cache[0][0],
@@ -144,18 +146,20 @@ def gamma_ns(order, mode, n, nf):
                 full_sx_cache[3][0],
             ]
         )
+        sx_new = full_sx_cache[-1]
     else:
         sx = harmonics.sx(n, max_weight=order[0] + 1)
+        sx_new = harmonics.compute_additional_sx_cache(n, False, False, sx)
     # 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 mode == 10101:
-            gamma_ns_1 = as2.gamma_nsp(n, nf, sx)
+            gamma_ns_1 = as2.gamma_nsp(n, nf, sx, sx_new)
         # To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here
         elif mode in [10201, 10200]:
-            gamma_ns_1 = as2.gamma_nsm(n, nf, sx)
+            gamma_ns_1 = as2.gamma_nsm(n, nf, sx, sx_new)
         else:
             raise NotImplementedError("Non-singlet sector is not implemented")
         gamma_ns[1] = gamma_ns_1
@@ -208,7 +212,9 @@ def gamma_singlet(order, n, nf):
     """
     # cache the s-es
     if order[0] >= 4:
-        full_sx_cache = harmonics.compute_cache(n, 5, is_singlet=False)
+        full_sx_cache = harmonics.compute_cache(
+            n, 5, is_singlet=False, additional=True, qed=False
+        )
         sx = np.array(
             [
                 full_sx_cache[0][0],
@@ -217,16 +223,18 @@ def gamma_singlet(order, n, nf):
                 full_sx_cache[3][0],
             ]
         )
+        sx_new = full_sx_cache[-1]
     elif order[0] >= 3:
         # here we need only S1,S2,S3,S4
         sx = harmonics.sx(n, max_weight=order[0] + 1)
+        sx_new = harmonics.compute_additional_sx_cache(n, True, False, sx)
     else:
         sx = harmonics.sx(n, max_weight=order[0])
-
+        sx_new = harmonics.compute_additional_sx_cache(n, True, False, sx)
     gamma_s = np.zeros((order[0], 2, 2), np.complex_)
     gamma_s[0] = as1.gamma_singlet(n, sx[0], nf)
     if order[0] >= 2:
-        gamma_s[1] = as2.gamma_singlet(n, nf, sx)
+        gamma_s[1] = as2.gamma_singlet(n, nf, sx, sx_new)
     if order[0] >= 3:
         gamma_s[2] = as3.gamma_singlet(n, nf, sx)
     if order[0] >= 4:
@@ -279,21 +287,21 @@ def gamma_ns_qed(order, mode, n, nf):
         sx = harmonics.sx(n, max_weight=max_weight + 1)
     else:
         sx = harmonics.sx(n, max_weight=3)
-    sx_ns_qed = harmonics.compute_qed_ns_cache(n, sx[0])
+    sx_new = harmonics.compute_additional_sx_cache(n, False, True, sx)
     # now combine
     gamma_ns = np.zeros((order[0] + 1, order[1] + 1), np.complex_)
     gamma_ns[1, 0] = as1.gamma_ns(n, sx[0])
     gamma_ns[0, 1] = choose_ns_ad_aem1(mode, n, sx)
-    gamma_ns[1, 1] = choose_ns_ad_as1aem1(mode, n, sx, sx_ns_qed)
+    gamma_ns[1, 1] = choose_ns_ad_as1aem1(mode, n, sx, sx_new)
     # NLO and beyond
     if order[0] >= 2:
         if mode in [10102, 10103]:
-            gamma_ns[2, 0] = as2.gamma_nsp(n, nf, sx)
+            gamma_ns[2, 0] = as2.gamma_nsp(n, nf, sx, sx_new)
         # To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here
         else:
-            gamma_ns[2, 0] = as2.gamma_nsm(n, nf, sx)
+            gamma_ns[2, 0] = as2.gamma_nsm(n, nf, sx, sx_new)
     if order[1] >= 2:
-        gamma_ns[0, 2] = choose_ns_ad_aem2(mode, n, nf, sx, sx_ns_qed)
+        gamma_ns[0, 2] = choose_ns_ad_aem2(mode, n, nf, sx, sx_new)
     # NNLO and beyond
     if order[0] >= 3:
         if mode in [10102, 10103]:
@@ -331,7 +339,7 @@ def choose_ns_ad_aem1(mode, n, sx):
 
 
 @nb.njit(cache=True)
-def choose_ns_ad_as1aem1(mode, n, sx, sx_ns_qed):
+def choose_ns_ad_as1aem1(mode, n, sx, sx_new):
     r"""
     Select the non-singlet anomalous dimension at O(as1aem1) with the correct charge factor.
 
@@ -350,17 +358,17 @@ def choose_ns_ad_as1aem1(mode, n, sx, sx_ns_qed):
             non-singlet anomalous dimensions
     """
     if mode == 10102:
-        return constants.eu2 * as1aem1.gamma_nsp(n, sx, sx_ns_qed)
+        return constants.eu2 * as1aem1.gamma_nsp(n, sx, sx_new)
     elif mode == 10103:
-        return constants.ed2 * as1aem1.gamma_nsp(n, sx, sx_ns_qed)
+        return constants.ed2 * as1aem1.gamma_nsp(n, sx, sx_new)
     elif mode == 10202:
-        return constants.eu2 * as1aem1.gamma_nsm(n, sx, sx_ns_qed)
+        return constants.eu2 * as1aem1.gamma_nsm(n, sx, sx_new)
     elif mode == 10203:
-        return constants.ed2 * as1aem1.gamma_nsm(n, sx, sx_ns_qed)
+        return constants.ed2 * as1aem1.gamma_nsm(n, sx, sx_new)
 
 
 @nb.njit(cache=True)
-def choose_ns_ad_aem2(mode, n, nf, sx, sx_ns_qed):
+def choose_ns_ad_aem2(mode, n, nf, sx, sx_new):
     r"""
     Select the non-singlet anomalous dimension at O(aem2) with the correct charge factor.
 
@@ -379,13 +387,13 @@ def choose_ns_ad_aem2(mode, n, nf, sx, sx_ns_qed):
             non-singlet anomalous dimensions
     """
     if mode == 10102:
-        return constants.eu2 * aem2.gamma_nspu(n, nf, sx, sx_ns_qed)
+        return constants.eu2 * aem2.gamma_nspu(n, nf, sx, sx_new)
     elif mode == 10103:
-        return constants.ed2 * aem2.gamma_nspd(n, nf, sx, sx_ns_qed)
+        return constants.ed2 * aem2.gamma_nspd(n, nf, sx, sx_new)
     elif mode == 10202:
-        return constants.eu2 * aem2.gamma_nsmu(n, nf, sx, sx_ns_qed)
+        return constants.eu2 * aem2.gamma_nsmu(n, nf, sx, sx_new)
     elif mode == 10203:
-        return constants.ed2 * aem2.gamma_nsmd(n, nf, sx, sx_ns_qed)
+        return constants.ed2 * aem2.gamma_nsmd(n, nf, sx, sx_new)
 
 
 @nb.njit(cache=True)
@@ -423,15 +431,15 @@ def gamma_singlet_qed(order, n, nf):
         sx = harmonics.sx(n, max_weight=max_weight + 1)
     else:
         sx = harmonics.sx(n, max_weight=3)
-    sx_ns_qed = harmonics.compute_qed_ns_cache(n, sx[0])
+    sx_new = harmonics.compute_additional_sx_cache(n, False, True, sx)
     gamma_s = np.zeros((order[0] + 1, order[1] + 1, 4, 4), np.complex_)
     gamma_s[1, 0] = as1.gamma_singlet_qed(n, sx[0], nf)
     gamma_s[0, 1] = aem1.gamma_singlet(n, nf, sx)
-    gamma_s[1, 1] = as1aem1.gamma_singlet(n, nf, sx, sx_ns_qed)
+    gamma_s[1, 1] = as1aem1.gamma_singlet(n, nf, sx, sx_new)
     if order[0] >= 2:
-        gamma_s[2, 0] = as2.gamma_singlet_qed(n, nf, sx)
+        gamma_s[2, 0] = as2.gamma_singlet_qed(n, nf, sx, sx_new)
     if order[1] >= 2:
-        gamma_s[0, 2] = aem2.gamma_singlet(n, nf, sx, sx_ns_qed)
+        gamma_s[0, 2] = aem2.gamma_singlet(n, nf, sx, sx_new)
     if order[0] >= 3:
         gamma_s[3, 0] = as3.gamma_singlet_qed(n, nf, sx)
     return gamma_s
@@ -472,15 +480,15 @@ def gamma_valence_qed(order, n, nf):
         sx = harmonics.sx(n, max_weight=max_weight + 1)
     else:
         sx = harmonics.sx(n, max_weight=3)
-    sx_ns_qed = harmonics.compute_qed_ns_cache(n, sx[0])
+    sx_new = harmonics.compute_additional_sx_cache(n, False, True, sx)
     gamma_v = np.zeros((order[0] + 1, order[1] + 1, 2, 2), np.complex_)
     gamma_v[1, 0] = as1.gamma_valence_qed(n, sx[0])
     gamma_v[0, 1] = aem1.gamma_valence(n, nf, sx)
-    gamma_v[1, 1] = as1aem1.gamma_valence(n, nf, sx, sx_ns_qed)
+    gamma_v[1, 1] = as1aem1.gamma_valence(n, nf, sx, sx_new)
     if order[0] >= 2:
-        gamma_v[2, 0] = as2.gamma_valence_qed(n, nf, sx)
+        gamma_v[2, 0] = as2.gamma_valence_qed(n, nf, sx, sx_new)
     if order[1] >= 2:
-        gamma_v[0, 2] = aem2.gamma_valence(n, nf, sx, sx_ns_qed)
+        gamma_v[0, 2] = aem2.gamma_valence(n, nf, sx, sx_new)
     if order[0] >= 3:
         gamma_v[3, 0] = as3.gamma_valence_qed(n, nf, sx)
     return gamma_v
diff --git a/src/eko/anomalous_dimensions/aem2.py b/src/eko/anomalous_dimensions/aem2.py
index 05bb033ad..77b9dd44e 100644
--- a/src/eko/anomalous_dimensions/aem2.py
+++ b/src/eko/anomalous_dimensions/aem2.py
@@ -151,7 +151,7 @@ def gamma_phd(N, nf, sx):
 
 
 @nb.njit(cache=True)
-def gamma_nspu(N, nf, sx, sx_ns_qed):
+def gamma_nspu(N, nf, sx, sx_new):
     r"""Compute the :math:`O(a_{em}^2)` singlet-like non-singlet anomalous dimension for up quarks.
 
     Implements sum of Eqs. (57-58) of :cite:`deFlorian:2016gvk` for q=u.
@@ -184,13 +184,11 @@ def gamma_nspu(N, nf, sx, sx_ns_qed):
         - 80.0 / 9.0 * S1
         + 16.0 / 3.0 * S2
     ) * eSigma2
-    return (
-        constants.eu2 * as1aem1.gamma_nsp(N, sx, sx_ns_qed) / constants.CF / 2.0 + tmp
-    )
+    return constants.eu2 * as1aem1.gamma_nsp(N, sx, sx_new) / constants.CF / 2.0 + tmp
 
 
 @nb.njit(cache=True)
-def gamma_nspd(N, nf, sx, sx_ns_qed):
+def gamma_nspd(N, nf, sx, sx_new):
     r"""Compute the :math:`O(a_{em}^2)` singlet-like non-singlet anomalous dimension for down quarks.
 
     Implements sum of Eqs. (57-58) of :cite:`deFlorian:2016gvk` for q=d.
@@ -223,13 +221,11 @@ def gamma_nspd(N, nf, sx, sx_ns_qed):
         - 80.0 / 9.0 * S1
         + 16.0 / 3.0 * S2
     ) * eSigma2
-    return (
-        constants.ed2 * as1aem1.gamma_nsp(N, sx, sx_ns_qed) / constants.CF / 2.0 + tmp
-    )
+    return constants.ed2 * as1aem1.gamma_nsp(N, sx, sx_new) / constants.CF / 2.0 + tmp
 
 
 @nb.njit(cache=True)
-def gamma_nsmu(N, nf, sx, sx_ns_qed):
+def gamma_nsmu(N, nf, sx, sx_new):
     r"""Compute the :math:`O(a_{em}^2)` valence-like non-singlet anomalous dimension for up quarks.
 
     Implements difference between Eqs. (57-58) of :cite:`deFlorian:2016gvk` for q=u.
@@ -262,13 +258,11 @@ def gamma_nsmu(N, nf, sx, sx_ns_qed):
         - 80.0 / 9.0 * S1
         + 16.0 / 3.0 * S2
     ) * eSigma2
-    return (
-        constants.eu2 * as1aem1.gamma_nsm(N, sx, sx_ns_qed) / constants.CF / 2.0 + tmp
-    )
+    return constants.eu2 * as1aem1.gamma_nsm(N, sx, sx_new) / constants.CF / 2.0 + tmp
 
 
 @nb.njit(cache=True)
-def gamma_nsmd(N, nf, sx, sx_ns_qed):
+def gamma_nsmd(N, nf, sx, sx_new):
     r"""Compute the :math:`O(a_{em}^2)` valence-like non-singlet anomalous dimension for down quarks.
 
     Implements difference between Eqs. (57-58) of :cite:`deFlorian:2016gvk` for q=d.
@@ -301,9 +295,7 @@ def gamma_nsmd(N, nf, sx, sx_ns_qed):
         - 80.0 / 9.0 * S1
         + 16.0 / 3.0 * S2
     ) * eSigma2
-    return (
-        constants.ed2 * as1aem1.gamma_nsm(N, sx, sx_ns_qed) / constants.CF / 2.0 + tmp
-    )
+    return constants.ed2 * as1aem1.gamma_nsm(N, sx, sx_new) / constants.CF / 2.0 + tmp
 
 
 @nb.njit(cache=True)
@@ -336,7 +328,7 @@ def gamma_ps(N, nf):
 
 
 @nb.njit(cache=True)
-def gamma_singlet(N, nf, sx, sx_ns_qed):
+def gamma_singlet(N, nf, sx, sx_new):
     r"""Compute the :math:`O(a_{em}^2)` singlet sector.
 
     Parameters
@@ -363,8 +355,8 @@ def gamma_singlet(N, nf, sx, sx_ns_qed):
     gamma_ph_d = gamma_phd(N, nf, sx)
     gamma_u_ph = gamma_uph(N, nf, sx)
     gamma_d_ph = gamma_dph(N, nf, sx)
-    gamma_ns_p_u = gamma_nspu(N, nf, sx, sx_ns_qed)
-    gamma_ns_p_d = gamma_nspd(N, nf, sx, sx_ns_qed)
+    gamma_ns_p_u = gamma_nspu(N, nf, sx, sx_new)
+    gamma_ns_p_d = gamma_nspd(N, nf, sx, sx_new)
     gamma_pure_singlet = gamma_ps(N, nf)
     gamma_S_02 = np.array(
         [
@@ -408,7 +400,7 @@ def gamma_singlet(N, nf, sx, sx_ns_qed):
 
 
 @nb.njit(cache=True)
-def gamma_valence(N, nf, sx, sx_ns_qed):
+def gamma_valence(N, nf, sx, sx_new):
     r"""Compute the :math:`O(a_{em}^2)` valence sector.
 
     Parameters
@@ -429,8 +421,8 @@ def gamma_valence(N, nf, sx, sx_ns_qed):
     nd = nf - nu
     vu = nu / nf
     vd = nd / nf
-    gamma_ns_m_u = gamma_nsmu(N, nf, sx, sx_ns_qed)
-    gamma_ns_m_d = gamma_nsmd(N, nf, sx, sx_ns_qed)
+    gamma_ns_m_u = gamma_nsmu(N, nf, sx, sx_new)
+    gamma_ns_m_d = gamma_nsmd(N, nf, sx, sx_new)
     gamma_V_02 = np.array(
         [
             [
diff --git a/src/eko/anomalous_dimensions/as1aem1.py b/src/eko/anomalous_dimensions/as1aem1.py
index f3f9b8ba6..12aaf30d5 100644
--- a/src/eko/anomalous_dimensions/as1aem1.py
+++ b/src/eko/anomalous_dimensions/as1aem1.py
@@ -241,7 +241,7 @@ def gamma_gg():
 
 
 @nb.njit(cache=True)
-def gamma_nsp(N, sx, sx_ns_qed):
+def gamma_nsp(N, sx, sx_new):
     r"""Compute the :math:`O(a_s^1a_{em}^1)` singlet-like non-singlet anomalous dimension.
 
     Implements sum of Eqs. (33-34) of :cite:`deFlorian:2015ujt`.
@@ -263,14 +263,14 @@ def gamma_nsp(N, sx, sx_ns_qed):
     S1 = sx[0]
     S2 = sx[1]
     S3 = sx[2]
-    S1h = sx_ns_qed[0]
-    S2h = sx_ns_qed[1]
-    S3h = sx_ns_qed[2]
-    S1p1h = sx_ns_qed[3]
-    S2p1h = sx_ns_qed[4]
-    S3p1h = sx_ns_qed[5]
-    g3N = sx_ns_qed[6]
-    g3Np2 = sx_ns_qed[7]
+    S1h = sx_new[0]
+    S2h = sx_new[1]
+    S3h = sx_new[2]
+    g3N = sx_new[3]
+    g3Np2 = sx_new[7]
+    S1p1h = sx_new[8]
+    S2p1h = sx_new[9]
+    S3p1h = sx_new[10]
     result = (
         +32.0 * zeta2 * S1h
         - 32.0 * zeta2 * S1p1h
@@ -308,7 +308,7 @@ def gamma_nsp(N, sx, sx_ns_qed):
 
 
 @nb.njit(cache=True)
-def gamma_nsm(N, sx, sx_ns_qed):
+def gamma_nsm(N, sx, sx_new):
     r"""Compute the :math:`O(a_s^1a_{em}^1)` valence-like non-singlet anomalous dimension.
 
     Implements difference between Eqs. (33-34) of :cite:`deFlorian:2015ujt`.
@@ -330,14 +330,14 @@ def gamma_nsm(N, sx, sx_ns_qed):
     S1 = sx[0]
     S2 = sx[1]
     S3 = sx[2]
-    S1h = sx_ns_qed[0]
-    S2h = sx_ns_qed[1]
-    S3h = sx_ns_qed[2]
-    S1p1h = sx_ns_qed[3]
-    S2p1h = sx_ns_qed[4]
-    S3p1h = sx_ns_qed[5]
-    g3N = sx_ns_qed[6]
-    g3Np2 = sx_ns_qed[7]
+    S1h = sx_new[0]
+    S2h = sx_new[1]
+    S3h = sx_new[2]
+    g3N = sx_new[3]
+    g3Np2 = sx_new[7]
+    S1p1h = sx_new[8]
+    S2p1h = sx_new[9]
+    S3p1h = sx_new[10]
     result = (
         -32.0 * zeta2 * S1h
         - 8.0 / (N + N**2) * S2h
@@ -371,7 +371,7 @@ def gamma_nsm(N, sx, sx_ns_qed):
 
 
 @nb.njit(cache=True)
-def gamma_singlet(N, nf, sx, sx_ns_qed):
+def gamma_singlet(N, nf, sx, sx_new):
     r"""Compute the :math:`O(a_s^1a_{em}^1)` singlet sector.
 
     Parameters
@@ -397,7 +397,7 @@ def gamma_singlet(N, nf, sx, sx_ns_qed):
     gamma_ph_q = gamma_phq(N, sx)
     gamma_q_g = gamma_qg(N, nf, sx)
     gamma_q_ph = gamma_qph(N, nf, sx)
-    gamma_ns_p = gamma_nsp(N, sx, sx_ns_qed)
+    gamma_ns_p = gamma_nsp(N, sx, sx_new)
     gamma_S_11 = np.array(
         [
             [
@@ -431,7 +431,7 @@ def gamma_singlet(N, nf, sx, sx_ns_qed):
 
 
 @nb.njit(cache=True)
-def gamma_valence(N, nf, sx, sx_ns_qed):
+def gamma_valence(N, nf, sx, sx_new):
     r"""Compute the :math:`O(a_s^1a_{em}^1)` valence sector.
 
     Parameters
@@ -459,4 +459,4 @@ def gamma_valence(N, nf, sx, sx_ns_qed):
         ],
         np.complex_,
     )
-    return gamma_V_11 * gamma_nsm(N, sx, sx_ns_qed)
+    return gamma_V_11 * gamma_nsm(N, sx, sx_new)
diff --git a/src/eko/anomalous_dimensions/as2.py b/src/eko/anomalous_dimensions/as2.py
index 9d263eefa..2792e288f 100644
--- a/src/eko/anomalous_dimensions/as2.py
+++ b/src/eko/anomalous_dimensions/as2.py
@@ -13,7 +13,7 @@
 
 
 @nb.njit(cache=True)
-def gamma_nsm(n, nf, sx):
+def gamma_nsm(n, nf, sx, sx_new):
     r"""Compute the |NLO| valence-like non-singlet anomalous dimension.
 
     Implements Eq. (3.5) of :cite:`Moch:2004pa`.
@@ -37,12 +37,10 @@ def gamma_nsm(n, nf, sx):
     S2 = sx[1]
     # Here, Sp refers to S' ("s-prime") (german: "s-strich" or in Pegasus language: SSTR)
     # of :cite:`Gluck:1989ze` and NOT to the Spence function a.k.a. dilogarithm
-    # TODO : these harmonic sums are computed also for the QED sector then we can use
-    # the ones that are passed to the O(as1aem1) anomalous dimensions
-    Sp1m = harmonics.S1((n - 1) / 2)
-    Sp2m = harmonics.S2((n - 1) / 2)
-    Sp3m = harmonics.S3((n - 1) / 2)
-    g3n = harmonics.g_functions.mellin_g3(n, S1)
+    g3n = sx_new[3]
+    Sp1m = sx_new[4]
+    Sp2m = sx_new[5]
+    Sp3m = sx_new[6]
     # fmt: off
     gqq1m_cfca = 16*g3n - (144 + n*(1 + n)*(156 + n*(340 + n*(655 + 51*n*(2 + n)))))/(18.*np.power(n,3)*np.power(1 + n,3)) + (-14.666666666666666 + 8/n - 8/(1 + n))*S2 - (4*Sp2m)/(n + np.power(n,2)) + S1*(29.77777777777778 + 16/np.power(n,2) - 16*S2 + 8*Sp2m) + 2*Sp3m + 10*zeta3 + zeta2*(16*S1 - 16*Sp1m - (16*(1 + n*log2))/n) # pylint: disable=line-too-long
     gqq1m_cfcf = -32*g3n + (24 - n*(-32 + 3*n*(-8 + n*(3 + n)*(3 + np.power(n,2)))))/(2.*np.power(n,3)*np.power(1 + n,3)) + (12 - 8/n + 8/(1 + n))*S2 + S1*(-24/np.power(n,2) - 8/np.power(1 + n,2) + 16*S2 - 16*Sp2m) + (8*Sp2m)/(n + np.power(n,2)) - 4*Sp3m - 20*zeta3 + zeta2*(-32*S1 + 32*Sp1m + 32*(1/n + log2)) # pylint: disable=line-too-long
@@ -57,7 +55,7 @@ def gamma_nsm(n, nf, sx):
 
 
 @nb.njit(cache=True)
-def gamma_nsp(n, nf, sx):
+def gamma_nsp(n, nf, sx, sx_new):
     r"""Compute the |NLO| singlet-like non-singlet anomalous dimension.
 
     Implements Eq. (3.5) of :cite:`Moch:2004pa`.
@@ -79,10 +77,10 @@ def gamma_nsp(n, nf, sx):
     """
     S1 = sx[0]
     S2 = sx[1]
-    Sp1p = harmonics.S1(n / 2)
-    Sp2p = harmonics.S2(n / 2)
-    Sp3p = harmonics.S3(n / 2)
-    g3n = harmonics.g_functions.mellin_g3(n, S1)
+    Sp1p = sx_new[0]
+    Sp2p = sx_new[1]
+    Sp3p = sx_new[2]
+    g3n = sx_new[3]
     # fmt: off
     gqq1p_cfca = -16*g3n + (132 - n*(340 + n*(655 + 51*n*(2 + n))))/(18.*np.power(n,2)*np.power(1 + n,2)) + (-14.666666666666666 + 8/n - 8/(1 + n))*S2 - (4*Sp2p)/(n + np.power(n,2)) + S1*(29.77777777777778 - 16/np.power(n,2) - 16*S2 + 8*Sp2p) + 2*Sp3p + 10*zeta3 + zeta2*(16*S1 - 16*Sp1p + 16*(1/n - log2)) # pylint: disable=line-too-long
     gqq1p_cfcf = 32*g3n - (8 + n*(32 + n*(40 + 3*n*(3 + n)*(3 + np.power(n,2)))))/(2.*np.power(n,3)*np.power(1 + n,3)) + (12 - 8/n + 8/(1 + n))*S2 + S1*(40/np.power(n,2) - 8/np.power(1 + n,2) + 16*S2 - 16*Sp2p) + (8*Sp2p)/(n + np.power(n,2)) - 4*Sp3p - 20*zeta3 + zeta2*(-32*S1 + 32*Sp1p + 32*(-(1/n) + log2)) # pylint: disable=line-too-long
@@ -123,7 +121,7 @@ def gamma_ps(n, nf):
 
 
 @nb.njit(cache=True)
-def gamma_qg(n, nf, sx):
+def gamma_qg(n, nf, sx, sx_new):
     r"""Compute the |NLO| quark-gluon singlet anomalous dimension.
 
     Implements Eq. (3.7) of :cite:`Vogt:2004mw`.
@@ -145,7 +143,7 @@ def gamma_qg(n, nf, sx):
     """
     S1 = sx[0]
     S2 = sx[1]
-    Sp2p = harmonics.S2(n / 2)
+    Sp2p = sx_new[1]
     # fmt: off
     gqg1_nfca = (-4*(16 + n*(64 + n*(104 + n*(128 + n*(85 + n*(36 + n*(25 + n*(15 + n*(6 + n))))))))))/((-1 + n)*np.power(n,3)*np.power(1 + n,3)*np.power(2 + n,3)) - (16*(3 + 2*n)*S1)/np.power(2 + 3*n + np.power(n,2),2) + (4*(2 + n + np.power(n,2))*np.power(S1,2))/(n*(2 + 3*n + np.power(n,2))) - (4*(2 + n + np.power(n,2))*S2)/(n*(2 + 3*n + np.power(n,2))) + (4*(2 + n + np.power(n,2))*Sp2p)/(n*(2 + 3*n + np.power(n,2))) # pylint: disable=line-too-long
     gqg1_nfcf = (-2*(4 + n*(8 + n*(1 + n)*(25 + n*(26 + 5*n*(2 + n))))))/(np.power(n,3)*np.power(1 + n,3)*(2 + n)) + (8*S1)/np.power(n,2) - (4*(2 + n + np.power(n,2))*np.power(S1,2))/(n*(2 + 3*n + np.power(n,2))) + (4*(2 + n + np.power(n,2))*S2)/(n*(2 + 3*n + np.power(n,2))) # pylint: disable=line-too-long
@@ -157,7 +155,7 @@ def gamma_qg(n, nf, sx):
 
 
 @nb.njit(cache=True)
-def gamma_gq(n, nf, sx):
+def gamma_gq(n, nf, sx, sx_new):
     r"""Compute the |NLO| gluon-quark singlet anomalous dimension.
 
     Implements Eq. (3.8) of :cite:`Vogt:2004mw`.
@@ -179,7 +177,7 @@ def gamma_gq(n, nf, sx):
     """
     S1 = sx[0]
     S2 = sx[1]
-    Sp2p = harmonics.S2(n / 2)
+    Sp2p = sx_new[1]
     # fmt: off
     ggq1_cfcf = (-8 + 2*n*(-12 + n*(-1 + n*(28 + n*(43 + 6*n*(5 + 2*n))))))/((-1 + n)*np.power(n,3)*np.power(1 + n,3)) - (4*(10 + n*(17 + n*(8 + 5*n)))*S1)/((-1 + n)*n*np.power(1 + n,2)) + (4*(2 + n + np.power(n,2))*np.power(S1,2))/(n*(-1 + np.power(n,2))) + (4*(2 + n + np.power(n,2))*S2)/(n*(-1 + np.power(n,2))) # pylint: disable=line-too-long
     ggq1_cfca = (-4*(144 + n*(432 + n*(-152 + n*(-1304 + n*(-1031 + n*(695 + n*(1678 + n*(1400 + n*(621 + 109*n))))))))))/(9.*np.power(n,3)*np.power(1 + n,3)*np.power(-2 + n + np.power(n,2),2)) + (4*(-12 + n*(-22 + 41*n + 17*np.power(n,3)))*S1)/(3.*np.power(-1 + n,2)*np.power(n,2)*(1 + n)) + ((8 + 4*n + 4*np.power(n,2))*np.power(S1,2))/(n - np.power(n,3)) + ((8 + 4*n + 4*np.power(n,2))*S2)/(n - np.power(n,3)) + (4*(2 + n + np.power(n,2))*Sp2p)/(n*(-1 + np.power(n,2))) # pylint: disable=line-too-long
@@ -194,7 +192,7 @@ def gamma_gq(n, nf, sx):
 
 
 @nb.njit(cache=True)
-def gamma_gg(n, nf, sx):
+def gamma_gg(n, nf, sx, sx_new):
     r"""Compute the |NLO| gluon-gluon singlet anomalous dimension.
 
     Implements Eq. (3.9) of :cite:`Vogt:2004mw`.
@@ -215,10 +213,10 @@ def gamma_gg(n, nf, sx):
         :math:`\\gamma_{gg}^{(1)}(N)`
     """
     S1 = sx[0]
-    Sp1p = harmonics.S1(n / 2)
-    Sp2p = harmonics.S2(n / 2)
-    Sp3p = harmonics.S3(n / 2)
-    g3n = harmonics.g_functions.mellin_g3(n, S1)
+    Sp1p = sx_new[0]
+    Sp2p = sx_new[1]
+    Sp3p = sx_new[2]
+    g3n = sx_new[3]
     # fmt: off
     ggg1_caca = 16*g3n - (2*(576 + n*(1488 + n*(560 + n*(-1248 + n*(-1384 + n*(1663 + n*(4514 + n*(4744 + n*(3030 + n*(1225 + 48*n*(7 + n))))))))))))/(9.*np.power(-1 + n,2)*np.power(n,3)*np.power(1 + n,3)*np.power(2 + n,3)) + S1*(29.77777777777778 + 16/np.power(-1 + n,2) + 16/np.power(1 + n,2) - 16/np.power(2 + n,2) - 8*Sp2p) + (16*(1 + n + np.power(n,2))*Sp2p)/(n*(1 + n)*(-2 + n + np.power(n,2))) - 2*Sp3p - 10*zeta3 + zeta2*(-16*S1 + 16*Sp1p + 16*(-(1/n) + log2)) # pylint: disable=line-too-long
     ggg1_canf = (8*(6 + n*(1 + n)*(28 + n*(1 + n)*(13 + 3*n*(1 + n)))))/(9.*np.power(n,2)*np.power(1 + n,2)*(-2 + n + np.power(n,2))) - (40*S1)/9. # pylint: disable=line-too-long
@@ -231,7 +229,7 @@ def gamma_gg(n, nf, sx):
 
 
 @nb.njit(cache=True)
-def gamma_singlet(n, nf, sx):
+def gamma_singlet(n, nf, sx, sx_new):
     r"""Compute the next-leading-order singlet anomalous dimension matrix.
 
     .. math::
@@ -262,16 +260,19 @@ def gamma_singlet(n, nf, sx):
     gamma_gq : :math:`\\gamma_{gq}^{(1)}`
     gamma_gg : :math:`\\gamma_{gg}^{(1)}`
     """
-    gamma_qq = gamma_nsp(n, nf, sx) + gamma_ps(n, nf)
+    gamma_qq = gamma_nsp(n, nf, sx, sx_new) + gamma_ps(n, nf)
     gamma_S_0 = np.array(
-        [[gamma_qq, gamma_qg(n, nf, sx)], [gamma_gq(n, nf, sx), gamma_gg(n, nf, sx)]],
+        [
+            [gamma_qq, gamma_qg(n, nf, sx, sx_new)],
+            [gamma_gq(n, nf, sx, sx_new), gamma_gg(n, nf, sx, sx_new)],
+        ],
         np.complex_,
     )
     return gamma_S_0
 
 
 @nb.njit(cache=True)
-def gamma_singlet_qed(N, nf, sx):
+def gamma_singlet_qed(N, nf, sx, sx_new):
     r"""Compute the leading-order singlet anomalous dimension matrix.
 
     .. math::
@@ -303,13 +304,18 @@ def gamma_singlet_qed(N, nf, sx):
     gamma_gq : :math:`\\gamma_{gq}^{(0)}`
     gamma_gg : :math:`\\gamma_{gg}^{(0)}`
     """
-    gamma_ns_p = gamma_nsp(N, nf, sx)
+    gamma_ns_p = gamma_nsp(N, nf, sx, sx_new)
     gamma_qq = gamma_ns_p + gamma_ps(N, nf)
     gamma_S = np.array(
         [
-            [gamma_gg(N, nf, sx), 0.0 + 0.0j, gamma_gq(N, nf, sx), 0.0 + 0.0j],
+            [
+                gamma_gg(N, nf, sx, sx_new),
+                0.0 + 0.0j,
+                gamma_gq(N, nf, sx, sx_new),
+                0.0 + 0.0j,
+            ],
             [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
-            [gamma_qg(N, nf, sx), 0.0 + 0.0j, gamma_qq, 0.0 + 0.0j],
+            [gamma_qg(N, nf, sx, sx_new), 0.0 + 0.0j, gamma_qq, 0.0 + 0.0j],
             [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, gamma_ns_p],
         ],
         np.complex_,
@@ -318,7 +324,7 @@ def gamma_singlet_qed(N, nf, sx):
 
 
 @nb.njit(cache=True)
-def gamma_valence_qed(N, nf, sx):
+def gamma_valence_qed(N, nf, sx, sx_new):
     r"""Compute the leading-order valence anomalous dimension matrix.
 
     .. math::
@@ -355,4 +361,4 @@ def gamma_valence_qed(N, nf, sx):
         ],
         np.complex_,
     )
-    return gamma_V * gamma_nsm(N, nf, sx)
+    return gamma_V * gamma_nsm(N, nf, sx, sx_new)
diff --git a/src/eko/harmonics/__init__.py b/src/eko/harmonics/__init__.py
index 02cbe78c2..5d07de63d 100644
--- a/src/eko/harmonics/__init__.py
+++ b/src/eko/harmonics/__init__.py
@@ -6,6 +6,7 @@
 import numpy as np
 
 from . import g_functions, polygamma
+from .constants import log2, zeta2, zeta3
 from .w1 import S1, Sm1
 from .w2 import S2, Sm2
 from .w3 import S3, S21, S2m1, Sm2m1, Sm3, Sm21
@@ -180,7 +181,7 @@ def s4x(n, sx, smx, is_singlet):
 
 
 @nb.njit(cache=True)
-def compute_cache(n, max_weight, is_singlet):
+def compute_cache(n, max_weight, is_singlet, additional=False, qed=False):
     r"""Get the harmonics sums cache.
 
     Parameters
@@ -217,19 +218,48 @@ def compute_cache(n, max_weight, is_singlet):
         sx[2, 1:-2] = s3x(n, sx[:, 0], sx[:, -1], is_singlet)
         sx[3, 1:-1] = s4x(n, sx[:, 0], sx[:, -1], is_singlet)
     # return list of list keeping the non zero values
-    return [[el for el in sx_list if el != 0] for sx_list in sx]
+    tmp = [[el for el in sx_list if el != 0] for sx_list in sx]
+    if additional:
+        tmp.append(
+            compute_additional_sx_cache(
+                n,
+                is_singlet,
+                qed,
+                sx[:, 0],
+                sx[:, -1],
+                sx[2, 1 if max_weight == 3 else 3],
+            )
+        )
+    return tmp
 
 
 @nb.njit(cache=True)
-def compute_qed_ns_cache(n, s1):
+def compute_additional_sx_cache(n, is_singlet, qed, sx, sminus=None, sm21=None):
     r"""Get the harmonics sums cache needed for the qed non singlet AD.
 
+    We have three cases :
+    - 1 : for qcd singlet we need only S_1(n/2), S_2(n/2), S_{3}(n/2) and g_3(n)
+    - 2 : for qcd non singlet we need also S_1((n-1)/2), S_2((n-1)/2), S_{3}((n-1)/2)
+    - 3 : for qed we need also g_3(n+2), S_1((n-1)/2), S_2((n-1)/2), S_{3}((n-1)/2)
+    These three cases are called with (qcd_singlet, qed) = :
+    - 1 : (True, False)
+    - 2 : (False, False)
+    - 3 : (False, True)
+
     Parameters
     ----------
     n : complex
         Mellin moment
-    s1 : float
-        harmonic sum :math:`S_1(N)`
+    qcd_singlet : bool
+        True if singlet for pure QCD evolution
+    qed : bool
+        True if QED evolution activated
+    s1 : list
+        harmonic sum :math:`S_i(N)`
+    s1 : list
+        harmonic sum :math:`S_{-i}(N)`
+    sm21 : complex
+        harmonic sum :math:`S_{-2,1}(N)`
 
     Returns
     -------
@@ -240,28 +270,90 @@ def compute_qed_ns_cache(n, s1):
             [S_1(n/2),
             S_2(n/2),
             S_{3}(n/2),
+            g_3(n),
+            S_1((n-1)/2),
+            S_2((n-1)/2),
+            S_{3}((n-1)/2),
+            g_3(n+2),
             S_1((n+1)/2),
             S_2((n+1)/2),
-            S_{3}((n+1)/2),
-            g_3(n),
-            g_3(n+2)]
+            S_{3}((n+1)/2)]
+
+    """
+    if sminus is None:
+        sm1 = None  # needed later
+        s1h = S1(n / 2.0)
+        s2h = S2(n / 2.0)
+        s3h = S3(n / 2.0)
+        g3N = g_functions.mellin_g3(n, sx[0])
+        sx_new = [s1h, s2h, s3h, g3N]
+        if not is_singlet:
+            s1_nm1_h = S1((n - 1.0) / 2.0)
+            s2_nm1_h = S2((n - 1.0) / 2.0)
+            s3_nm1_h = S3((n - 1.0) / 2.0)
+            sx_new.extend([s1_nm1_h, s2_nm1_h, s3_nm1_h])
+    else:
+        sm1 = sminus[0]  # needed later
+        if is_singlet:
+            s1h = sminus[0] + sx[0]
+            s2h = 2 * (sminus[1] + sx[1])
+            s3h = 4 * (sminus[2] + sx[2])
+            g3N = -(
+                zeta2 * (sminus[0] - 1 / n)
+                - 5.0 / 8 * zeta3
+                + zeta2 * log2
+                - (sm21 - 1 / n**2 * sx[0])
+            )
+            sx_new = [s1h, s2h, s3h, g3N]
+        else:
+            s1h = S1(n / 2.0)
+            s2h = S2(n / 2.0)
+            s3h = S3(n / 2.0)
+            g3N = (
+                zeta2 * (sminus[0] + 1 / n)
+                - 5.0 / 8 * zeta3
+                + zeta2 * log2
+                - (sm21 + 1 / n**2 * sx[0])
+            )
+            s1_nm1_h = sminus[0] + sx[0]
+            s2_nm1_h = 2 * (sminus[1] + sx[1])
+            s3_nm1_h = 4 * (sminus[2] + sx[2])
+            sx_new = [s1h, s2h, s3h, g3N, s1_nm1_h, s2_nm1_h, s3_nm1_h]
+    if qed:
+        g3Np2 = g3np2(n, is_singlet, sx[0], sm21, sm1)
+        s1_np1_h = polygamma.recursive_harmonic_sum(s1_nm1_h, (n - 1.0) / 2.0, 1, 1)
+        s2_np1_h = polygamma.recursive_harmonic_sum(s2_nm1_h, (n - 1.0) / 2.0, 1, 2)
+        s3_np1_h = polygamma.recursive_harmonic_sum(s3_nm1_h, (n - 1.0) / 2.0, 1, 3)
+        sx_new.extend([g3Np2, s1_np1_h, s2_np1_h, s3_np1_h])
+    return sx_new
+
+
+@nb.njit(cache=True)
+def g3np2(n, is_singlet, s1, sm21=None, sm1=None):
+    """
+    Compute g3(n+2).
 
+    Parameters
+    ----------
+    n : complex
+        Mellin moment
+    is_singlet : bool
+        True if singlet for pure QCD evolution
+    qed : bool
+        True if QED evolution activated
+    s1 : complex
+        harmonic sum :math:`S_1(N)`
+
+    Returns
+    -------
+    g3(n+2) : complex
     """
-    s1h = S1(n / 2.0)
-    sx_qed_ns = [s1h]
-    S2h = S2(n / 2.0)
-    sx_qed_ns.append(S2h)
-    S3h = S3(n / 2.0)
-    sx_qed_ns.append(S3h)
-    S1p1h = S1((n + 1.0) / 2.0)
-    sx_qed_ns.append(S1p1h)
-    S2p1h = S2((n + 1.0) / 2.0)
-    sx_qed_ns.append(S2p1h)
-    S3p1h = S3((n + 1.0) / 2.0)
-    sx_qed_ns.append(S3p1h)
-    g3N = g_functions.mellin_g3(n, s1)
-    sx_qed_ns.append(g3N)
-    S1p2 = polygamma.recursive_harmonic_sum(s1, n, 2, 1)
-    g3Np2 = g_functions.mellin_g3(n + 2.0, S1p2)
-    sx_qed_ns.append(g3Np2)
-    return sx_qed_ns
+    eta = 1 if is_singlet else -1
+    if sm21 is None:
+        s1_np2 = polygamma.recursive_harmonic_sum(s1, n, 2, 1)
+        return g_functions.mellin_g3(n + 2.0, s1_np2)
+    else:
+        s1_np1 = polygamma.recursive_harmonic_sum(s1, n, 1, 1)
+        sm1_np1 = sm1 - eta / (n + 1)
+        sm21_np1 = sm21 - eta / (n + 1) ** 2 * s1_np1
+        return -1 / eta * (zeta2 * sm1_np1 - 5.0 / 8 * zeta3 + zeta2 * log2 - sm21_np1)
diff --git a/tests/eko/anomalous_dimensions/test_aem2.py b/tests/eko/anomalous_dimensions/test_aem2.py
index 362b46306..92ebb0788 100644
--- a/tests/eko/anomalous_dimensions/test_aem2.py
+++ b/tests/eko/anomalous_dimensions/test_aem2.py
@@ -10,7 +10,7 @@ def test_number_conservation():
     # number
     N = complex(1.0, 0.0)
     sx = h.sx(N, 3)
-    sx_ns_qed = h.compute_qed_ns_cache(N, sx[0])
+    sx_ns_qed = h.compute_additional_sx_cache(N, False, True, sx)
     for NF in range(2, 6 + 1):
         np.testing.assert_almost_equal(
             ad.aem2.gamma_nsmu(N, NF, sx, sx_ns_qed), 0, decimal=4
@@ -39,7 +39,7 @@ def test_quark_momentum_conservation():
     # quark momentum
     N = complex(2.0, 0.0)
     sx = h.sx(N, 3)
-    sx_ns_qed = h.compute_qed_ns_cache(N, sx[0])
+    sx_ns_qed = h.compute_additional_sx_cache(N, False, True, sx)
     NF = 6
     NU = constants.uplike_flavors(NF)
     ND = NF - NU
diff --git a/tests/eko/anomalous_dimensions/test_as1aem1.py b/tests/eko/anomalous_dimensions/test_as1aem1.py
index c0db6dd6b..58cb11e05 100644
--- a/tests/eko/anomalous_dimensions/test_as1aem1.py
+++ b/tests/eko/anomalous_dimensions/test_as1aem1.py
@@ -11,7 +11,7 @@ def test_number_conservation():
     # number
     N = complex(1.0, 0.0)
     sx = h.sx(N, 3)
-    sx_ns_qed = h.compute_qed_ns_cache(N, sx[0])
+    sx_ns_qed = h.compute_additional_sx_cache(N, False, True, sx)
     np.testing.assert_almost_equal(ad.as1aem1.gamma_nsm(N, sx, sx_ns_qed), 0, decimal=4)
 
 
@@ -53,7 +53,7 @@ def test_quark_momentum_conservation():
     # quark momentum
     N = complex(2.0, 0.0)
     sx = h.sx(N, 3)
-    sx_ns_qed = h.compute_qed_ns_cache(N, sx[0])
+    sx_ns_qed = h.compute_additional_sx_cache(N, False, True, sx)
     np.testing.assert_almost_equal(
         ad.as1aem1.gamma_nsp(N, sx, sx_ns_qed)
         + ad.as1aem1.gamma_gq(N, sx)
diff --git a/tests/eko/anomalous_dimensions/test_as2.py b/tests/eko/anomalous_dimensions/test_as2.py
index 01d7f2f28..6562833ea 100644
--- a/tests/eko/anomalous_dimensions/test_as2.py
+++ b/tests/eko/anomalous_dimensions/test_as2.py
@@ -11,10 +11,12 @@
 def test_gamma_1():
     # number conservation
     sx_n1 = h.sx(1, 2)
-    np.testing.assert_allclose(ad_as2.gamma_nsm(1, NF, sx_n1), 0.0, atol=2e-6)
+    sx_new1 = h.compute_additional_sx_cache(1, False, False, sx_n1)
+    np.testing.assert_allclose(ad_as2.gamma_nsm(1, NF, sx_n1, sx_new1), 0.0, atol=2e-6)
 
     sx_n2 = h.sx(2, 2)
-    gS1 = ad_as2.gamma_singlet(2, NF, sx_n2)
+    sx_new2 = h.compute_additional_sx_cache(2, True, False, sx_n2)
+    gS1 = ad_as2.gamma_singlet(2, NF, sx_n2, sx_new2)
     # gluon momentum conservation
     # the CA*NF term seems to be tough to compute, so raise the constraint ...
     np.testing.assert_allclose(gS1[0, 1] + gS1[1, 1], 0.0, atol=4e-5)
@@ -26,7 +28,7 @@ def test_gamma_1():
     # reference values are obtained from MMa
     # non-singlet sector
     np.testing.assert_allclose(
-        ad_as2.gamma_nsp(2, NF, sx_n2),
+        ad_as2.gamma_nsp(2, NF, sx_n2, sx_new2),
         (-112.0 * const.CF + 376.0 * const.CA - 64.0 * NF) * const.CF / 27.0,
     )
     # singlet sector
@@ -39,8 +41,9 @@ def test_gamma_1():
     )  # gq
 
     # add additional point at (analytical) continuation point
+    sx_new2 = h.compute_additional_sx_cache(2, False, False, sx_n2)
     np.testing.assert_allclose(
-        ad_as2.gamma_nsm(2, NF, sx_n2),
+        ad_as2.gamma_nsm(2, NF, sx_n2, sx_new2),
         (
             (34.0 / 27.0 * (-47.0 + 6 * np.pi**2) - 16.0 * h.constants.zeta3)
             * const.CF
@@ -52,8 +55,10 @@ def test_gamma_1():
     )
     sx_n3 = h.sx(3, 2)
     sx_n4 = h.sx(4, 2)
+    sx_new3 = h.compute_additional_sx_cache(3, True, False, sx_n3)
+    sx_new4 = h.compute_additional_sx_cache(4, True, False, sx_n4)
     np.testing.assert_allclose(
-        ad_as2.gamma_nsp(3, NF, sx_n3),
+        ad_as2.gamma_nsp(3, NF, sx_n3, sx_new3),
         (
             (-34487.0 / 432.0 + 86.0 * np.pi**2 / 9.0 - 16.0 * h.constants.zeta3)
             * const.CF
@@ -64,7 +69,7 @@ def test_gamma_1():
         * const.CF,
     )
     np.testing.assert_allclose(ad_as2.gamma_ps(3, NF), -1391.0 * const.CF * NF / 5400.0)
-    gS1 = ad_as2.gamma_singlet(3, NF, sx_n3)
+    gS1 = ad_as2.gamma_singlet(3, NF, sx_n3, sx_new3)
     np.testing.assert_allclose(
         gS1[1, 0],
         (
@@ -84,7 +89,7 @@ def test_gamma_1():
         ),
         rtol=6e-7,
     )  # gg
-    gS1 = ad_as2.gamma_singlet(4, NF, sx_n4)
+    gS1 = ad_as2.gamma_singlet(4, NF, sx_n4, sx_new4)
     np.testing.assert_allclose(
         gS1[0, 1], (-56317.0 / 18000.0 * const.CF + 16387.0 / 9000.0 * const.CA) * NF
     )  # qg
@@ -92,7 +97,7 @@ def test_gamma_1():
     const.update_colors(4)
 
     np.testing.assert_allclose(const.CA, 4.0)
-    gS1 = ad_as2.gamma_singlet(3, NF, sx_n3)
+    gS1 = ad_as2.gamma_singlet(3, NF, sx_n3, sx_new3)
     np.testing.assert_allclose(
         gS1[1, 0],
         (
@@ -112,7 +117,7 @@ def test_gamma_1():
         ),
         rtol=6e-7,
     )  # gg
-    gS1 = ad_as2.gamma_singlet(4, NF, sx_n4)
+    gS1 = ad_as2.gamma_singlet(4, NF, sx_n4, sx_new4)
     np.testing.assert_allclose(
         gS1[0, 1], (-56317.0 / 18000.0 * const.CF + 16387.0 / 9000.0 * const.CA) * NF
     )  # qg
diff --git a/tests/eko/anomalous_dimensions/test_init.py b/tests/eko/anomalous_dimensions/test_init.py
index dce1066ff..b133ca9b7 100644
--- a/tests/eko/anomalous_dimensions/test_init.py
+++ b/tests/eko/anomalous_dimensions/test_init.py
@@ -195,11 +195,12 @@ def test_dim_singlet():
     nf = 3
     N = 2
     sx = harmonics.sx(N, max_weight=3 + 1)
+    sx_new = harmonics.compute_additional_sx_cache(N, False, True, sx)
     gamma_singlet = ad.gamma_singlet_qed((3, 2), N, nf)
     assert gamma_singlet.shape == (4, 3, 4, 4)
     gamma_singlet_as1 = as1.gamma_singlet_qed(N, sx[0], nf)
     assert gamma_singlet_as1.shape == (4, 4)
-    gamma_singlet_as2 = as2.gamma_singlet_qed(N, nf, sx)
+    gamma_singlet_as2 = as2.gamma_singlet_qed(N, nf, sx, sx_new)
     assert gamma_singlet_as2.shape == (4, 4)
     gamma_singlet_as3 = as3.gamma_singlet_qed(N, nf, sx)
     assert gamma_singlet_as3.shape == (4, 4)
@@ -209,11 +210,12 @@ def test_dim_valence():
     nf = 3
     N = 2
     sx = harmonics.sx(N, max_weight=3 + 1)
+    sx_new = harmonics.compute_additional_sx_cache(N, False, False, sx)
     gamma_valence = ad.gamma_valence_qed((3, 2), N, nf)
     assert gamma_valence.shape == (4, 3, 2, 2)
     gamma_valence_as1 = as1.gamma_valence_qed(N, sx[0])
     assert gamma_valence_as1.shape == (2, 2)
-    gamma_valence_as2 = as2.gamma_valence_qed(N, nf, sx)
+    gamma_valence_as2 = as2.gamma_valence_qed(N, nf, sx, sx_new)
     assert gamma_valence_as2.shape == (2, 2)
     gamma_valence_as3 = as3.gamma_valence_qed(N, nf, sx)
     assert gamma_valence_as3.shape == (2, 2)
diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py
index db9174c91..59f4fbacd 100644
--- a/tests/eko/evolution_operator/test_init.py
+++ b/tests/eko/evolution_operator/test_init.py
@@ -28,183 +28,40 @@ def test_quad_ker(monkeypatch):
     monkeypatch.setattr(ns, "dispatcher", lambda *args: 1.0)
     monkeypatch.setattr(qed_ns, "dispatcher", lambda *args: 1.0)
     monkeypatch.setattr(s, "dispatcher", lambda *args: np.identity(2))
-    for is_log in [True, False]:
-        res_ns = quad_ker(
-            u=0,
-            order=(1, 0),
-            mode0=br.non_singlet_pids_map["ns+"],
-            mode1=0,
-            method="",
-            is_log=is_log,
-            logx=0.0,
-            areas=np.zeros(3),
-            as1=1,
-            as0=2,
-            as_raw=1,
-            aem_list=[0.00058],
-            alphaem_running=False,
-            nf=3,
-            L=0,
-            ev_op_iterations=0,
-            ev_op_max_order=(0, 0),
-            sv_mode=1,
-            is_threshold=False,
-        )
-        np.testing.assert_allclose(res_ns, 0.0)
-        res_ns = quad_ker(
-            u=0,
-            order=(3, 1),
-            mode0=br.non_singlet_pids_map["ns+u"],
-            mode1=0,
-            method="",
-            is_log=is_log,
-            logx=0.0,
-            areas=np.zeros(3),
-            as1=1,
-            as0=2,
-            as_raw=1,
-            aem_list=[0.00058],
-            alphaem_running=False,
-            nf=3,
-            L=0,
-            ev_op_iterations=0,
-            ev_op_max_order=(0, 0),
-            sv_mode=1,
-            is_threshold=False,
-        )
-        np.testing.assert_allclose(res_ns, 0.0)
-        res_s = quad_ker(
-            u=0,
-            order=(1, 0),
-            mode0=100,
-            mode1=100,
-            method="",
-            is_log=is_log,
-            logx=0.123,
-            areas=np.zeros(3),
-            as1=1,
-            as0=2,
-            as_raw=1,
-            aem_list=[0.00058],
-            alphaem_running=False,
-            nf=3,
-            L=0,
-            ev_op_iterations=0,
-            ev_op_max_order=(0, 0),
-            sv_mode=1,
-            is_threshold=False,
-        )
-        np.testing.assert_allclose(res_s, 1.0)
-        res_s = quad_ker(
-            u=0,
-            order=(1, 1),
-            mode0=100,
-            mode1=100,
-            method="iterate-exact",
-            is_log=is_log,
-            logx=0.123,
-            areas=np.zeros(3),
-            as1=1,
-            as0=2,
-            as_raw=1,
-            aem_list=[0.00058],
-            alphaem_running=False,
-            nf=3,
-            L=0,
-            ev_op_iterations=0,
-            ev_op_max_order=(0, 0),
-            sv_mode=1,
-            is_threshold=False,
-        )
-        np.testing.assert_allclose(res_s, 1.0)
-        res_s = quad_ker(
-            u=0,
-            order=(1, 0),
-            mode0=100,
-            mode1=21,
-            method="",
-            is_log=is_log,
-            logx=0.0,
-            areas=np.zeros(3),
-            as1=1,
-            as0=2,
-            as_raw=1,
-            aem_list=[0.00058],
-            alphaem_running=False,
-            nf=3,
-            L=0,
-            ev_op_iterations=0,
-            ev_op_max_order=(0, 0),
-            sv_mode=1,
-            is_threshold=False,
-        )
-        np.testing.assert_allclose(res_s, 0.0)
-        res_s = quad_ker(
-            u=0,
-            order=(1, 1),
-            mode0=100,
-            mode1=21,
-            method="iterate-exact",
-            is_log=is_log,
-            logx=0.0,
-            areas=np.zeros(3),
-            as1=1,
-            as0=2,
-            as_raw=1,
-            aem_list=[0.00058],
-            alphaem_running=False,
-            nf=3,
-            L=0,
-            ev_op_iterations=0,
-            ev_op_max_order=(0, 0),
-            sv_mode=1,
-            is_threshold=False,
-        )
-        np.testing.assert_allclose(res_s, 0.0)
-        res_v = quad_ker(
-            u=0,
-            order=(1, 1),
-            mode0=10200,
-            mode1=10200,
-            method="iterate-exact",
-            is_log=is_log,
-            logx=0.123,
-            areas=np.zeros(3),
-            as1=1,
-            as0=2,
-            as_raw=1,
-            aem_list=[0.00058],
-            alphaem_running=False,
-            nf=3,
-            L=0,
-            ev_op_iterations=0,
-            ev_op_max_order=(0, 0),
-            sv_mode=1,
-            is_threshold=False,
-        )
-        np.testing.assert_allclose(res_v, 1.0)
-        res_v = quad_ker(
-            u=0,
-            order=(1, 1),
-            mode0=10200,
-            mode1=10204,
-            method="iterate-exact",
-            is_log=is_log,
-            logx=0.123,
-            areas=np.zeros(3),
-            as1=1,
-            as0=2,
-            as_raw=1,
-            aem_list=[0.00058],
-            alphaem_running=False,
-            nf=3,
-            L=0,
-            ev_op_iterations=0,
-            ev_op_max_order=(0, 0),
-            sv_mode=1,
-            is_threshold=False,
-        )
-        np.testing.assert_allclose(res_v, 0.0)
+    params = [
+        ((1, 0), br.non_singlet_pids_map["ns+"], 0, "", 0.0, 0.0),
+        ((3, 1), br.non_singlet_pids_map["ns+u"], 0, "", 0.0, 0.0),
+        ((1, 0), 100, 100, "", 0.123, 1.0),
+        ((1, 0), 100, 21, "", 0.0, 0.0),
+        ((1, 1), 100, 100, "iterate-exact", 0.123, 1.0),
+        ((1, 1), 100, 21, "iterate-exact", 0.123, 0.0),
+        ((1, 1), 10200, 10200, "iterate-exact", 0.123, 1.0),
+        ((1, 1), 10200, 10204, "iterate-exact", 0.123, 0.0),
+    ]
+    for p in params:
+        for is_log in [True, False]:
+            res_ns = quad_ker(
+                u=0,
+                order=p[0],
+                mode0=p[1],
+                mode1=p[2],
+                method=p[3],
+                is_log=is_log,
+                logx=p[4],
+                areas=np.zeros(3),
+                as1=1,
+                as0=2,
+                as_raw=1,
+                aem_list=[0.00058],
+                alphaem_running=False,
+                nf=3,
+                L=0,
+                ev_op_iterations=0,
+                ev_op_max_order=(0, 0),
+                sv_mode=1,
+                is_threshold=False,
+            )
+            np.testing.assert_allclose(res_ns, p[5])
     for label in [(br.non_singlet_pids_map["ns+"], 0), (100, 100)]:
         for sv in [2, 3]:
             res_sv = quad_ker(
@@ -287,14 +144,23 @@ def test_quad_ker(monkeypatch):
     np.testing.assert_allclose(res_ns, 0.0)
 
 
+class FakeCoupling:
+    def __init__(self):
+        self.alphaem_running = None
+        self.q2_ref = 0.0
+
+    def a(self, scale_to=None, fact_scale=None, nf_to=None):
+        return (0.1, 0.01)
+
+    def compute_aem_as(self, aem_ref, as_from, as_to, nf):
+        return aem_ref
+
+
+fake_managers = {"couplings": FakeCoupling()}
+
+
 class TestOperator:
-    def test_labels(self, theory_ffns, operator_card, tmp_path):
-        tcard: TheoryCard = theory_ffns(3)
-        tcard.xif = 2.0
-        ocard: OperatorCard = operator_card
-        ocard.configs.scvar_method = ScaleVariationsMethod.EXPONENTIATED
-        r = eko.runner.legacy.Runner(tcard, ocard, path=tmp_path / "eko.tar")
-        g = r.op_grid
+    def test_labels(self):
         o = Operator(
             dict(
                 order=(3, 0),
@@ -303,7 +169,7 @@ def test_labels(self, theory_ffns, operator_card, tmp_path):
                 n_integration_cores=1,
                 ModSV=None,
             ),
-            g.managers,
+            fake_managers,
             3,
             1,
             2,
@@ -317,20 +183,14 @@ def test_labels(self, theory_ffns, operator_card, tmp_path):
                 n_integration_cores=1,
                 ModSV=None,
             ),
-            g.managers,
+            fake_managers,
             3,
             1,
             2,
         )
         assert sorted(o.labels) == []
 
-    def test_labels_qed(self, theory_ffns, operator_card, tmp_path):
-        tcard: TheoryCard = theory_ffns(3)
-        tcard.xif = 2.0
-        ocard: OperatorCard = operator_card
-        ocard.configs.scvar_method = ScaleVariationsMethod.EXPONENTIATED
-        r = eko.runner.legacy.Runner(tcard, ocard, path=tmp_path / "eko.tar")
-        g = r.op_grid
+    def test_labels_qed(self):
         o = Operator(
             dict(
                 order=(3, 1),
@@ -340,7 +200,7 @@ def test_labels_qed(self, theory_ffns, operator_card, tmp_path):
                 ModSV=None,
                 ev_op_iterations=1,
             ),
-            g.managers,
+            fake_managers,
             3,
             1,
             2,
@@ -355,24 +215,18 @@ def test_labels_qed(self, theory_ffns, operator_card, tmp_path):
                 ModSV=None,
                 ev_op_iterations=1,
             ),
-            g.managers,
+            fake_managers,
             3,
             1,
             2,
         )
         assert sorted(o.labels) == []
 
-    def test_n_pools(self, theory_ffns, operator_card, tmp_path):
+    def test_n_pools(self):
         excluded_cores = 3
         # make sure we actually have more the those cores (e.g. on github we don't)
         if os.cpu_count() <= excluded_cores:
             return
-        tcard: TheoryCard = theory_ffns(3)
-        tcard.xif = 2.0
-        ocard: OperatorCard = operator_card
-        ocard.configs.scvar_method = ScaleVariationsMethod.EXPONENTIATED
-        r = eko.runner.legacy.Runner(tcard, ocard, path=tmp_path / "eko.tar")
-        g = r.op_grid
         o = Operator(
             dict(
                 order=(2, 0),
@@ -381,7 +235,7 @@ def test_n_pools(self, theory_ffns, operator_card, tmp_path):
                 n_integration_cores=-excluded_cores,
                 ModSV=None,
             ),
-            g.managers,
+            fake_managers,
             3,
             1,
             10,