diff --git a/src/eko/basis_rotation.py b/src/eko/basis_rotation.py
index 4b972592c..e31a7beb6 100644
--- a/src/eko/basis_rotation.py
+++ b/src/eko/basis_rotation.py
@@ -275,7 +275,7 @@
 }
 
 
-def ad_projector(ad_lab, nf, qed=False):
+def ad_projector(ad_lab, nf, qed):
     """
     Build a projector (as a numpy array) for the given anomalous dimension sector.
 
@@ -355,7 +355,7 @@ def select_light_flavors_uni_ev(ad_lab, nf):
         return map_ad_to_evolution[ad_lab]
 
 
-def ad_projectors(nf, qed=False):
+def ad_projectors(nf, qed):
     """
     Build projectors tensor (as a numpy array), collecting all the individual sector projectors.
 
diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index eab3c64de..29e67c19c 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -307,6 +307,7 @@ def quad_ker(
             ev_op_max_order,
             sv_mode,
             is_threshold,
+            n3lo_ad_variation,
         )
 
     # recombine everything
@@ -450,6 +451,7 @@ def quad_ker_qed(
     ev_op_max_order,
     sv_mode,
     is_threshold,
+    n3lo_ad_variation,
 ):
     """Raw evolution kernel inside quad.
 
@@ -489,6 +491,8 @@ def quad_ker_qed(
         scale variation mode, see `eko.scale_variations.Modes`
     is_threshold : boolean
         is this an itermediate threshold operator?
+    n3lo_ad_variation : tuple
+        |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)``
 
     Returns
     -------
@@ -497,7 +501,7 @@ def quad_ker_qed(
     """
     # compute the actual evolution kernel for QEDxQCD
     if ker_base.is_QEDsinglet:
-        gamma_s = ad_us.gamma_singlet_qed(order, ker_base.n, nf)
+        gamma_s = ad_us.gamma_singlet_qed(order, ker_base.n, nf, n3lo_ad_variation)
         # scale var exponentiated is directly applied on gamma
         if sv_mode == sv.Modes.exponentiated:
             gamma_s = sv.exponentiated.gamma_variation_qed(
@@ -622,7 +626,7 @@ def __init__(
         self.alphaem_running = self.managers["couplings"].alphaem_running
         if self.log_label == "Evolution":
             self.a = self.compute_a()
-            self.compute_aem_list()
+            self.as_list, self.a_half_list = self.compute_aem_list()
 
     @property
     def n_pools(self):
@@ -697,8 +701,8 @@ def compute_aem_list(self):
         """
         ev_op_iterations = self.config["ev_op_iterations"]
         if self.order[1] == 0:
-            self.as_list = np.array([self.a_s[0], self.a_s[1]])
-            self.a_half_list = np.zeros((ev_op_iterations, 2))
+            as_list = np.array([self.a_s[0], self.a_s[1]])
+            a_half = np.zeros((ev_op_iterations, 2))
         else:
             as0 = self.a_s[0]
             as1 = self.a_s[1]
@@ -718,7 +722,7 @@ def compute_aem_list(self):
             couplings = self.managers["couplings"]
             mu2_steps = utils.geomspace(self.q2_from, self.q2_to, 1 + ev_op_iterations)
             mu2_l = mu2_steps[0]
-            self.as_list = np.array(
+            as_list = np.array(
                 [
                     couplings.compute(
                         a_ref=a_start, nf=self.nf, scale_from=mu2_start, scale_to=mu2
@@ -734,7 +738,7 @@ def compute_aem_list(self):
                 )
                 a_half[step] = [a_s, aem]
                 mu2_l = mu2_h
-            self.a_half_list = a_half
+        return as_list, a_half
 
     @property
     def labels(self):
diff --git a/src/eko/evolution_operator/flavors.py b/src/eko/evolution_operator/flavors.py
index 7ad7b5801..fac94a608 100644
--- a/src/eko/evolution_operator/flavors.py
+++ b/src/eko/evolution_operator/flavors.py
@@ -52,7 +52,7 @@ def pids_from_intrinsic_evol(label, nf, normalize):
     return weights
 
 
-def get_range(evol_labels, qed=False):
+def get_range(evol_labels, qed):
     """Determine the number of light and heavy flavors participating in the input and output.
 
     Here, we assume that the T distributions (e.g. T15) appears *always*
@@ -73,7 +73,7 @@ def get_range(evol_labels, qed=False):
     nf_in = 3
     nf_out = 3
 
-    def update(label, qed=False):
+    def update(label, qed):
         nf = 3
         if label[0] == "T":
             if not qed:
@@ -129,7 +129,7 @@ def rotate_pm_to_flavor(label):
     return l
 
 
-def rotate_matching(nf, qed=False, inverse=False):
+def rotate_matching(nf, qed, inverse=False):
     """Rotation between matching basis (with e.g. S,g,...V8 and c+,c-) and new true evolution basis (with S,g,...V8,T15,V15).
 
     Parameters
@@ -206,7 +206,7 @@ def rotate_matching(nf, qed=False, inverse=False):
     return l
 
 
-def rotate_matching_inverse(nf, qed=False):
+def rotate_matching_inverse(nf, qed):
     """Inverse rotation between matching basis (with e.g. S,g,...V8 and c+,c-) and new true evolution basis (with S,g,...V8,T15,V15).
 
     Parameters
diff --git a/src/eko/evolution_operator/matching_condition.py b/src/eko/evolution_operator/matching_condition.py
index 059617fbf..8df4775ff 100644
--- a/src/eko/evolution_operator/matching_condition.py
+++ b/src/eko/evolution_operator/matching_condition.py
@@ -21,7 +21,7 @@ def split_ad_to_evol_map(
         nf,
         q2_thr,
         intrinsic_range,
-        qed=False,
+        qed,
     ):
         """
         Create the instance from the |OME|.
diff --git a/src/eko/evolution_operator/physical.py b/src/eko/evolution_operator/physical.py
index c963e2a26..1ee501a2a 100644
--- a/src/eko/evolution_operator/physical.py
+++ b/src/eko/evolution_operator/physical.py
@@ -22,7 +22,7 @@ class PhysicalOperator(member.OperatorBase):
     """
 
     @classmethod
-    def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range, qed=False):
+    def ad_to_evol_map(cls, op_members, nf, q2_final, intrinsic_range, qed):
         """
         Obtain map between the 3-dimensional anomalous dimension basis and the 4-dimensional evolution basis.
 
diff --git a/src/eko/kernels/non_singlet.py b/src/eko/kernels/non_singlet.py
index 24706a794..23b2ab47a 100644
--- a/src/eko/kernels/non_singlet.py
+++ b/src/eko/kernels/non_singlet.py
@@ -186,7 +186,7 @@ def nnlo_expanded(gamma_ns, a1, a0, beta):
 
 
 @nb.njit(cache=True)
-def n3lo_expanded(gamma_ns, a1, a0, nf):
+def n3lo_expanded(gamma_ns, a1, a0, beta):
     """|N3LO| non-singlet expanded EKO.
 
     Parameters
@@ -206,12 +206,8 @@ def n3lo_expanded(gamma_ns, a1, a0, nf):
         |N3LO| non-singlet expanded EKO
 
     """
-    beta0 = beta.beta_qcd((2, 0), nf)
-    b_list = [
-        beta.b_qcd((3, 0), nf),
-        beta.b_qcd((4, 0), nf),
-        beta.b_qcd((5, 0), nf),
-    ]
+    beta0 = beta[0]
+    b_list = [betas / beta0 for betas in beta[1:]]
     j12 = ei.j12(a1, a0, beta0)
     j13 = as4_ei.j13_expanded(a1, a0, beta0, b_list)
     j23 = as4_ei.j23_expanded(a1, a0, beta0, b_list)
@@ -225,7 +221,7 @@ def n3lo_expanded(gamma_ns, a1, a0, nf):
 
 
 @nb.njit(cache=True)
-def n3lo_exact(gamma_ns, a1, a0, nf):
+def n3lo_exact(gamma_ns, a1, a0, beta):
     """|N3LO| non-singlet exact EKO.
 
     Parameters
@@ -245,12 +241,8 @@ def n3lo_exact(gamma_ns, a1, a0, nf):
         |N3LO| non-singlet exact EKO
 
     """
-    beta0 = beta.beta_qcd((2, 0), nf)
-    b_list = [
-        beta.b_qcd((3, 0), nf),
-        beta.b_qcd((4, 0), nf),
-        beta.b_qcd((5, 0), nf),
-    ]
+    beta0 = beta[0]
+    b_list = [betas / beta0 for betas in beta[1:]]
     roots = as4_ei.roots(b_list)
     j12 = ei.j12(a1, a0, beta0)
     j13 = as4_ei.j13_exact(a1, a0, beta0, b_list, roots)
@@ -420,6 +412,6 @@ def dispatcher(
             "decompose-expanded",
             "perturbative-expanded",
         ]:
-            return n3lo_expanded(gamma_ns, a1, a0, nf)
-        return n3lo_exact(gamma_ns, a1, a0, nf)
+            return n3lo_expanded(gamma_ns, a1, a0, betalist)
+        return n3lo_exact(gamma_ns, a1, a0, betalist)
     raise NotImplementedError("Selected order is not implemented")
diff --git a/src/eko/kernels/non_singlet_qed.py b/src/eko/kernels/non_singlet_qed.py
index 86ec1a9ff..36d4d3c35 100644
--- a/src/eko/kernels/non_singlet_qed.py
+++ b/src/eko/kernels/non_singlet_qed.py
@@ -120,6 +120,29 @@ def as3_exact(gamma_ns, a1, a0, beta):
     return ns.nnlo_exact(gamma_ns, a1, a0, beta)
 
 
+@nb.njit(cache=True)
+def as4_exact(gamma_ns, a1, a0, beta):
+    """O(as3aem1) non-singlet exact EKO.
+
+    Parameters
+    ----------
+    gamma_ns : numpy.ndarray
+        non-singlet anomalous dimensions
+    a1 : float
+        target coupling value
+    a0 : float
+        initial coupling value
+    beta : list
+        list of the values of the beta functions
+
+    Returns
+    -------
+    e_ns^3 : complex
+        O(as4aem1) non-singlet exact EKO
+    """
+    return ns.n3lo_exact(gamma_ns, a1, a0, beta)
+
+
 @nb.njit(cache=True)
 def dispatcher(
     order,
@@ -216,6 +239,8 @@ def fixed_alphaem_exact(order, gamma_ns, a1, a0, aem, nf, mu2_from, mu2_to):
         qcd_only = as2_exact(gamma_ns_list[1:], a1, a0, betalist)
     elif order[0] == 3:
         qcd_only = as3_exact(gamma_ns_list[1:], a1, a0, betalist)
+    elif order[0] == 4:
+        qcd_only = as4_exact(gamma_ns_list[1:], a1, a0, betalist)
     else:
         raise NotImplementedError("Selected order is not implemented")
     return qcd_only * apply_qed(gamma_ns_list[0], mu2_from, mu2_to)
diff --git a/src/eko/member.py b/src/eko/member.py
index f03b0e9b4..edf8f89a3 100644
--- a/src/eko/member.py
+++ b/src/eko/member.py
@@ -266,7 +266,7 @@ def operator_multiply(left, right, operation):
                     new_oms[new_key] += operation(l_op, r_op)
         return new_oms
 
-    def to_flavor_basis_tensor(self, qed: bool = False):
+    def to_flavor_basis_tensor(self, qed: bool):
         """Convert the computations into an rank 4 tensor.
 
         A sparse tensor defined with dot-notation (e.g. ``S.g``) is converted
diff --git a/src/ekobox/apply.py b/src/ekobox/apply.py
index 354461c0a..336900798 100644
--- a/src/ekobox/apply.py
+++ b/src/ekobox/apply.py
@@ -12,7 +12,6 @@ def apply_pdf(
     lhapdf_like,
     targetgrid=None,
     rotate_to_evolution_basis=False,
-    qed=False,
 ):
     """
     Apply all available operators to the input PDFs.
@@ -34,6 +33,7 @@ def apply_pdf(
         out_grid : dict
             output PDFs and their associated errors for the computed mu2grid
     """
+    qed = eko.theory_card.order[1] > 0
     if rotate_to_evolution_basis:
         if not qed:
             rotate_flavor_to_evolution = br.rotate_flavor_to_evolution
@@ -99,6 +99,7 @@ def apply_pdf_flavor(
         if error_final is not None:
             out_grid[ep]["errors"] = dict(zip(eko.bases.targetpids, error_final))
 
+    qed = eko.theory_card.order[1] > 0
     # rotate to evolution basis
     if flavor_rotation is not None:
         for q2, op in out_grid.items():
diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
index 6a3c6cf06..6358a5af8 100644
--- a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py
@@ -155,6 +155,11 @@ def gamma_ns_qed(order, mode, n, nf):
             gamma_ns[3, 0] = as3.gamma_nsp(n, nf, cache)
         elif mode in [10202, 10203]:
             gamma_ns[3, 0] = as3.gamma_nsm(n, nf, cache)
+    if order[0] >= 4:
+        if mode in [10102, 10103]:
+            gamma_ns[4, 0] = as4.gamma_nsp(n, nf, cache)
+        elif mode in [10202, 10203]:
+            gamma_ns[4, 0] = as4.gamma_nsm(n, nf, cache)
     return gamma_ns
 
 
@@ -250,7 +255,7 @@ def choose_ns_ad_aem2(mode, n, nf, cache):
 
 
 @nb.njit(cache=True)
-def gamma_singlet_qed(order, n, nf):
+def gamma_singlet_qed(order, n, nf, n3lo_ad_variation):
     r"""
     Compute the grid of the QED singlet anomalous dimensions matrices.
 
@@ -262,6 +267,8 @@ def gamma_singlet_qed(order, n, nf):
         Mellin variable
     nf : int
         Number of active flavors
+    n3lo_ad_variation : tuple
+        |N3LO| anomalous dimension variation ``(gg_var, gq_var, qg_var, qq_var)``
 
     Returns
     -------
@@ -280,6 +287,8 @@ def gamma_singlet_qed(order, n, nf):
         gamma_s[0, 2] = aem2.gamma_singlet(n, nf, cache)
     if order[0] >= 3:
         gamma_s[3, 0] = as3.gamma_singlet_qed(n, nf, cache)
+    if order[0] >= 4:
+        gamma_s[4, 0] = as4.gamma_singlet_qed(n, nf, cache, n3lo_ad_variation)
     return gamma_s
 
 
@@ -314,4 +323,6 @@ def gamma_valence_qed(order, n, nf):
         gamma_v[0, 2] = aem2.gamma_valence(n, nf, cache)
     if order[0] >= 3:
         gamma_v[3, 0] = as3.gamma_valence_qed(n, nf, cache)
+    if order[0] >= 4:
+        gamma_v[4, 0] = as4.gamma_valence_qed(n, nf, cache)
     return gamma_v
diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py
index 10c4cb57c..56c2c16ff 100644
--- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py
+++ b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/__init__.py
@@ -55,3 +55,84 @@ def gamma_singlet(N, nf, cache, variation):
         np.complex_,
     )
     return gamma_S_0
+
+
+@nb.njit(cache=True)
+def gamma_singlet_qed(N, nf, cache, variation):
+    r"""Compute the leading-order singlet anomalous dimension matrix for the unified evolution basis.
+
+    .. math::
+        \\gamma_S^{(3,0)} = \\left(\begin{array}{cccc}
+        \\gamma_{gg}^{(3,0)} & 0 & \\gamma_{gq}^{(3,0)} & 0\\
+        0 & 0 & 0 & 0 \\
+        \\gamma_{qg}^{(3,0)} & 0 & \\gamma_{qq}^{(3,0)} & 0 \\
+        0 & 0 & 0 & \\gamma_{qq}^{(3,0)} \\
+        \\end{array}\right)
+
+    Parameters
+    ----------
+    N : complex
+        Mellin moment
+    nf : int
+        Number of active flavors
+    cache: numpy.ndarray
+        Harmonic sum cache
+
+    Returns
+    -------
+    numpy.ndarray
+        Leading-order singlet anomalous dimension matrix :math:`\\gamma_{S}^{(3,0)}(N)`
+
+    """
+    gamma_np_p = gamma_nsp(N, nf, cache)
+    gamma_qq = gamma_np_p + gamma_ps(N, nf, cache, variation[3])
+    gamma_S = np.array(
+        [
+            [
+                gamma_gg(N, nf, cache, variation[0]),
+                0.0 + 0.0j,
+                gamma_gq(N, nf, cache, variation[1]),
+                0.0 + 0.0j,
+            ],
+            [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j],
+            [gamma_qg(N, nf, cache, variation[2]), 0.0 + 0.0j, gamma_qq, 0.0 + 0.0j],
+            [0.0 + 0.0j, 0.0 + 0.0j, 0.0 + 0.0j, gamma_np_p],
+        ],
+        np.complex_,
+    )
+    return gamma_S
+
+
+@nb.njit(cache=True)
+def gamma_valence_qed(N, nf, cache):
+    r"""Compute the leading-order valence anomalous dimension matrix for the unified evolution basis.
+
+    .. math::
+        \\gamma_V^{(3,0)} = \\left(\begin{array}{cc}
+        \\gamma_{nsV}^{(3,0)} & 0\\
+        0 & \\gamma_{ns-}^{(3,0)}
+        \\end{array}\right)
+
+    Parameters
+    ----------
+    N : complex
+        Mellin moment
+    nf : int
+        Number of active flavors
+    cache: numpy.ndarray
+        Harmonic sum cache
+
+    Returns
+    -------
+    numpy.ndarray
+        Leading-order singlet anomalous dimension matrix :math:`\\gamma_{V}^{(3,0)}(N)`
+
+    """
+    gamma_V = np.array(
+        [
+            [gamma_nsv(N, nf, cache), 0.0],
+            [0.0, gamma_nsm(N, nf, cache)],
+        ],
+        np.complex_,
+    )
+    return gamma_V
diff --git a/tests/eko/evolution_operator/test_flavors.py b/tests/eko/evolution_operator/test_flavors.py
index 81bb35fc4..2f3761404 100644
--- a/tests/eko/evolution_operator/test_flavors.py
+++ b/tests/eko/evolution_operator/test_flavors.py
@@ -61,13 +61,15 @@ def get(d):
 
 
 def test_get_range():
-    assert (3, 3) == flavors.get_range([])
-    assert (3, 3) == flavors.get_range([member.MemberName(n) for n in ["S.S", "V3.V3"]])
+    assert (3, 3) == flavors.get_range([], False)
+    assert (3, 3) == flavors.get_range(
+        [member.MemberName(n) for n in ["S.S", "V3.V3"]], False
+    )
     assert (3, 4) == flavors.get_range(
-        [member.MemberName(n) for n in ["S.S", "V3.V3", "T15.S"]]
+        [member.MemberName(n) for n in ["S.S", "V3.V3", "T15.S"]], False
     )
     assert (4, 4) == flavors.get_range(
-        [member.MemberName(n) for n in ["S.S", "V3.V3", "T15.T15"]]
+        [member.MemberName(n) for n in ["S.S", "V3.V3", "T15.T15"]], False
     )
     assert (3, 3) == flavors.get_range(
         [member.MemberName(n) for n in ["S.S", "Td3.Td3"]], True
@@ -113,7 +115,7 @@ def test_rotate_pm_to_flavor():
 
 
 def test_rotate_matching():
-    m = flavors.rotate_matching(4)
+    m = flavors.rotate_matching(4, False)
     assert len(list(filter(lambda e: "c+" in e, m.keys()))) == 2
     assert len(list(filter(lambda e: "b-" in e, m.keys()))) == 1
 
@@ -156,8 +158,8 @@ def load(m):
         return mm
 
     for nf in range(4, 6 + 1):
-        m = load(flavors.rotate_matching(nf))
-        minv = load(flavors.rotate_matching_inverse(nf))
+        m = load(flavors.rotate_matching(nf, False))
+        minv = load(flavors.rotate_matching_inverse(nf, False))
         np.testing.assert_allclose(m @ minv, np.eye(len(br.evol_basis)), atol=1e-10)
 
 
diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py
index 642de5bb7..47d5e7001 100644
--- a/tests/eko/evolution_operator/test_init.py
+++ b/tests/eko/evolution_operator/test_init.py
@@ -382,7 +382,7 @@ def test_compute(self, monkeypatch, theory_ffns, operator_card, tmp_path):
                     err_msg=k,
                 )
 
-        for n in range(1, 3 + 1):
+        for n in range(1, 4 + 1):
             for qed in range(1, 2 + 1):
                 g.config["order"] = (n, qed)
                 o1 = Operator(g.config, g.managers, Segment(2.0, 2.0, 3))
diff --git a/tests/eko/evolution_operator/test_matching_condition.py b/tests/eko/evolution_operator/test_matching_condition.py
index 8cbd6b537..c854fb90c 100644
--- a/tests/eko/evolution_operator/test_matching_condition.py
+++ b/tests/eko/evolution_operator/test_matching_condition.py
@@ -38,7 +38,7 @@ def update_intrinsic_OME(self, ome):
 
     def test_split_ad_to_evol_map(self):
         ome = self.mkOME()
-        a = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [])
+        a = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [], False)
         triv_keys = [
             "V.V",
             "T3.T3",
@@ -64,7 +64,7 @@ def test_split_ad_to_evol_map(self):
             ome[(200, 200)].value,
         )
         # # if alpha is zero, nothing non-trivial should happen
-        b = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [])
+        b = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [], False)
         assert sorted(str(k) for k in b.op_members.keys()) == sorted(
             [*triv_keys, *keys3]
         )
@@ -74,7 +74,7 @@ def test_split_ad_to_evol_map(self):
         # )
         # nf=3 + IC
         self.update_intrinsic_OME(ome)
-        c = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [4])
+        c = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [4], False)
         assert sorted(str(k) for k in c.op_members.keys()) == sorted(
             [*triv_keys, *keys3, "S.c+", "g.c+", "c+.c+", "c-.c-"]
         )
@@ -83,7 +83,7 @@ def test_split_ad_to_evol_map(self):
             b.op_members[member.MemberName("V.V")].value,
         )
         # nf=3 + IB
-        d = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [5])
+        d = MatchingCondition.split_ad_to_evol_map(ome, 3, 1, [5], False)
         assert sorted(str(k) for k in d.op_members.keys()) == sorted(
             [*triv_keys, *keys3, "b+.b+", "b-.b-"]
         )
@@ -92,7 +92,7 @@ def test_split_ad_to_evol_map(self):
             np.eye(self.shape[0]),
         )
         # nf=4 + IB
-        d = MatchingCondition.split_ad_to_evol_map(ome, 4, 1, [5])
+        d = MatchingCondition.split_ad_to_evol_map(ome, 4, 1, [5], False)
         assert sorted(str(k) for k in d.op_members.keys()) == sorted(
             [
                 *triv_keys,
diff --git a/tests/eko/evolution_operator/test_physical.py b/tests/eko/evolution_operator/test_physical.py
index b4a9936b6..e9882de5b 100644
--- a/tests/eko/evolution_operator/test_physical.py
+++ b/tests/eko/evolution_operator/test_physical.py
@@ -138,7 +138,7 @@ def test_to_flavor_basis_tensor_ss(self):
             ),
             1,
         )
-        vt, _ = a.to_flavor_basis_tensor()
+        vt, _ = a.to_flavor_basis_tensor(False)
         np.testing.assert_allclose(vt[6, :, 6, :], vt[6, :, 5, :])
         np.testing.assert_allclose(vt[6, :, 6, :], vt[5, :, 6, :])
         np.testing.assert_allclose(vt[6, :, 6, :], SS.value[:, :] / (2 * 3))
@@ -158,7 +158,7 @@ def test_to_flavor_basis_tensor_gg(self):
             ),
             1,
         )
-        vt, _ = a.to_flavor_basis_tensor()
+        vt, _ = a.to_flavor_basis_tensor(False)
         np.testing.assert_allclose(vt[6, :, 6, :], 0)
         np.testing.assert_allclose(vt[7, :, 7, :], gg.value[:, :])
         np.testing.assert_allclose(vt[1, :, :, :], 0)
@@ -199,7 +199,7 @@ def test_to_flavor_basis_tensor_gg_qed(self):
             ),
             1,
         )
-        vt, _ = a.to_flavor_basis_tensor()
+        vt, _ = a.to_flavor_basis_tensor(False)
         np.testing.assert_allclose(vt[6, :, 6, :], 0)
         np.testing.assert_allclose(vt[7, :, 7, :], gg.value[:, :])
         np.testing.assert_allclose(vt[1, :, :, :], 0)
diff --git a/tests/eko/kernels/test_kernels_QEDns.py b/tests/eko/kernels/test_kernels_QEDns.py
index e84956084..e6179ba99 100644
--- a/tests/eko/kernels/test_kernels_QEDns.py
+++ b/tests/eko/kernels/test_kernels_QEDns.py
@@ -25,7 +25,7 @@ def test_zero():
     nf = 3
     ev_op_iterations = 2
     running_alpha = [True, False]
-    for qcd in range(1, 3 + 1):
+    for qcd in range(1, 4 + 1):
         for qed in range(1, 2 + 1):
             order = (qcd, qed)
             gamma_ns = (
@@ -73,7 +73,7 @@ def test_zero_true_gamma():
     for mode in br.non_singlet_pids_map.values():
         if mode in [10201, 10101, 10200]:
             continue
-        for qcd in range(1, 3 + 1):
+        for qcd in range(1, 4 + 1):
             for qed in range(1, 2 + 1):
                 order = (qcd, qed)
                 n = np.random.rand()
@@ -136,7 +136,8 @@ def test_ode():
     for mode in br.non_singlet_pids_map.values():
         if mode in [10201, 10101, 10200]:
             continue
-        for qcd in range(1, 3 + 1):
+        max_qcd = 4
+        for qcd in range(1, max_qcd + 1):
             for qed in range(1, 2 + 1):
                 order = (qcd, qed)
                 sc = Couplings(
@@ -154,14 +155,17 @@ def test_ode():
                     )
                     a1 = sc.a_s(mu2_to)
                     gamma_ns = (
-                        np.random.rand(3 + 1, 2 + 1) + np.random.rand(3 + 1, 2 + 1) * 1j
+                        np.random.rand(max_qcd + 1, 2 + 1)
+                        + np.random.rand(max_qcd + 1, 2 + 1) * 1j
                     )
                     gamma_ns[0, 0] = 0.0
                     gamma_ns[2, 1] = 0.0
                     gamma_ns[3, 1] = 0.0
+                    gamma_ns[4, 1] = 0.0
                     gamma_ns[1, 2] = 0.0
                     gamma_ns[2, 2] = 0.0
                     gamma_ns[3, 2] = 0.0
+                    gamma_ns[4, 2] = 0.0
                     gammatot = 0.0
                     for i in range(0, order[0] + 1):
                         for j in range(0, order[1] + 1):
@@ -226,7 +230,7 @@ def test_ode_true_gamma():
     for mode in br.non_singlet_pids_map.values():
         if mode in [10201, 10101, 10200]:
             continue
-        for qcd in range(1, 3 + 1):
+        for qcd in range(1, 4 + 1):
             for qed in range(1, 2 + 1):
                 order = (qcd, qed)
                 sc = Couplings(
@@ -302,9 +306,9 @@ def test_ode_true_gamma():
 
 def test_error():
     for running in [True, False]:
-        with pytest.raises(NotImplementedError):
+        with pytest.raises(ValueError):
             ns.dispatcher(
-                (4, 2),
+                (5, 2),
                 "iterate-exact",
                 np.random.rand(4, 3),
                 [0.1, 0.2],
diff --git a/tests/eko/kernels/test_kernels_QEDsinglet.py b/tests/eko/kernels/test_kernels_QEDsinglet.py
index 302bb941b..b953c053c 100644
--- a/tests/eko/kernels/test_kernels_QEDsinglet.py
+++ b/tests/eko/kernels/test_kernels_QEDsinglet.py
@@ -21,7 +21,7 @@ def test_zero(monkeypatch):
     nf = 3
     ev_op_iterations = 2
     ev_op_max_order = (3, 0)
-    for qcd in range(1, 3 + 1):
+    for qcd in range(1, 4 + 1):
         for qed in range(1, 2 + 1):
             order = (qcd, qed)
             gamma_s = (
@@ -75,11 +75,11 @@ def test_zero_true_gamma(monkeypatch):
     nf = 3
     ev_op_iterations = 2
     ev_op_max_order = (3, 0)
-    for qcd in range(1, 3 + 1):
+    for qcd in range(1, 4 + 1):
         for qed in range(1, 2 + 1):
             order = (qcd, qed)
             n = np.random.rand()
-            gamma_s = ad.gamma_singlet_qed(order, n, nf)
+            gamma_s = ad.gamma_singlet_qed(order, n, nf, (0, 0, 0, 0))
             # monkeypatch.setattr(
             #     ad,
             #     "exp_matrix_2D",
diff --git a/tests/eko/kernels/test_kernels_QEDvalence.py b/tests/eko/kernels/test_kernels_QEDvalence.py
index 95d0cbc4c..a27bc6ef6 100644
--- a/tests/eko/kernels/test_kernels_QEDvalence.py
+++ b/tests/eko/kernels/test_kernels_QEDvalence.py
@@ -21,7 +21,7 @@ def test_zero(monkeypatch):
     nf = 3
     ev_op_iterations = 2
     ev_op_max_order = (3, 0)
-    for qcd in range(1, 3 + 1):
+    for qcd in range(1, 4 + 1):
         for qed in range(1, 2 + 1):
             order = (qcd, qed)
             gamma_v = (
@@ -62,7 +62,7 @@ def test_zero_true_gamma(monkeypatch):
     nf = 3
     ev_op_iterations = 2
     ev_op_max_order = (3, 0)
-    for qcd in range(1, 3 + 1):
+    for qcd in range(1, 4 + 1):
         for qed in range(1, 2 + 1):
             order = (qcd, qed)
             n = np.random.rand()
diff --git a/tests/eko/test_basis_rotation.py b/tests/eko/test_basis_rotation.py
index 680ef065d..194669d46 100644
--- a/tests/eko/test_basis_rotation.py
+++ b/tests/eko/test_basis_rotation.py
@@ -9,19 +9,19 @@ def test_ad_projector():
     g = br.rotate_flavor_to_evolution[2]
     v3 = br.rotate_flavor_to_evolution[br.evol_basis.index("V3")]
 
-    s_to_s = br.ad_projector((100, 100), nf=6)
+    s_to_s = br.ad_projector((100, 100), nf=6, qed=False)
 
     np.testing.assert_allclose(s @ s_to_s, s)
     np.testing.assert_allclose(g @ s_to_s, 0.0)
     np.testing.assert_allclose(v3 @ s_to_s, 0.0)
 
-    g_to_s = br.ad_projector((21, 100), nf=6)
+    g_to_s = br.ad_projector((21, 100), nf=6, qed=False)
 
     np.testing.assert_allclose(s @ g_to_s, 0.0)
     np.testing.assert_allclose(g @ g_to_s, s)
     np.testing.assert_allclose(v3 @ g_to_s, 0.0)
 
-    ns_m = br.ad_projector((br.non_singlet_pids_map["ns-"], 0), nf=6)
+    ns_m = br.ad_projector((br.non_singlet_pids_map["ns-"], 0), nf=6, qed=False)
 
     np.testing.assert_allclose(s @ ns_m, 0.0, atol=1e-15)
     np.testing.assert_allclose(g @ ns_m, 0.0)
@@ -74,7 +74,7 @@ def test_ad_projectors():
     for nf in range(3, 6 + 1):
         diag = np.array([0] * (1 + 6 - nf) + [1] * (1 + 2 * nf) + [0] * (6 - nf))
         identity = np.diag(diag)
-        projs = br.ad_projectors(nf)
+        projs = br.ad_projectors(nf, qed=False)
 
         # sum over diagonal projectors form an identity
         np.testing.assert_allclose(
diff --git a/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py
index 2175eddcd..c7c955960 100644
--- a/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py
+++ b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py
@@ -208,7 +208,7 @@ def test_dim_singlet():
     nf = 3
     N = 2
     cache = h.cache.reset()
-    gamma_singlet = ad_us.gamma_singlet_qed((3, 2), N, nf)
+    gamma_singlet = ad_us.gamma_singlet_qed((3, 2), N, nf, (0, 0, 0, 0))
     assert gamma_singlet.shape == (4, 3, 4, 4)
     gamma_singlet_as1 = ad_us.as1.gamma_singlet_qed(N, cache, nf)
     assert gamma_singlet_as1.shape == (4, 4)