From eb0004ce3b57ee52766612fb89b48e5465a6e6a4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Tue, 9 Jan 2024 16:29:45 +0100
Subject: [PATCH 01/22] Pass nl as an argument in qed beta functions

---
 src/eko/beta.py      | 16 +++++++---------
 src/eko/constants.py |  2 --
 2 files changed, 7 insertions(+), 11 deletions(-)

diff --git a/src/eko/beta.py b/src/eko/beta.py
index 8ce6658f1..6d8b77ec2 100644
--- a/src/eko/beta.py
+++ b/src/eko/beta.py
@@ -32,7 +32,7 @@ def beta_qcd_as2(nf):
 
 
 @nb.njit(cache=True)
-def beta_qed_aem2(nf):
+def beta_qed_aem2(nf, nl=3):
     r"""Compute the first coefficient of the QED beta function.
 
     Implements Eq. (7) of :cite:`Surguladze:1996hx`.
@@ -50,7 +50,6 @@ def beta_qed_aem2(nf):
     """
     nu = constants.uplike_flavors(nf)
     nd = nf - nu
-    nl = 3  # TODO : pass nl as an argument??
     beta_qed_aem2 = (
         -4.0 / 3 * (nl + constants.NC * (nu * constants.eu2 + nd * constants.ed2))
     )
@@ -83,7 +82,7 @@ def beta_qcd_as3(nf):
 
 
 @nb.njit(cache=True)
-def beta_qed_aem3(nf):
+def beta_qed_aem3(nf, nl=3):
     r"""Compute the second coefficient of the QED beta function.
 
     Implements Eq. (7) of :cite:`Surguladze:1996hx`.
@@ -101,7 +100,6 @@ def beta_qed_aem3(nf):
     """
     nu = constants.uplike_flavors(nf)
     nd = nf - nu
-    nl = 3  # TODO : pass nl as an argument??
     beta_qed_aem3 = -4.0 * (
         nl + constants.NC * (nu * constants.eu2**2 + nd * constants.ed2**2)
     )
@@ -246,7 +244,7 @@ def beta_qcd(k, nf):
 
 
 @nb.njit(cache=True)
-def beta_qed(k, nf):
+def beta_qed(k, nf, nl=3):
     r"""Compute value of a beta_qed coefficients.
 
     Parameters
@@ -264,9 +262,9 @@ def beta_qed(k, nf):
     """
     beta_ = 0
     if k == (0, 2):
-        beta_ = beta_qed_aem2(nf)
+        beta_ = beta_qed_aem2(nf, nl)
     elif k == (0, 3):
-        beta_ = beta_qed_aem3(nf)
+        beta_ = beta_qed_aem3(nf, nl)
     elif k == (1, 2):
         beta_ = beta_qed_aem2as1(nf)
     else:
@@ -295,7 +293,7 @@ def b_qcd(k, nf):
 
 
 @nb.njit(cache=True)
-def b_qed(k, nf):
+def b_qed(k, nf, nl=3):
     r"""Compute b_qed coefficient.
 
     Parameters
@@ -311,4 +309,4 @@ def b_qed(k, nf):
             b_qed_k(nf)
 
     """
-    return beta_qed(k, nf) / beta_qed((0, 2), nf)
+    return beta_qed(k, nf, nl) / beta_qed((0, 2), nf, nl)
diff --git a/src/eko/constants.py b/src/eko/constants.py
index fbd0df459..c01197405 100644
--- a/src/eko/constants.py
+++ b/src/eko/constants.py
@@ -84,8 +84,6 @@ def uplike_flavors(nf):
     nu : int
 
     """
-    if nf not in range(2, 6 + 1):
-        raise NotImplementedError("Selected nf is not implemented")
     nu = nf // 2
     return nu
 

From ecec6f18e7cc14feca203a545a9d573334aa28d5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Tue, 9 Jan 2024 16:43:51 +0100
Subject: [PATCH 02/22] Restore check

---
 src/eko/constants.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/eko/constants.py b/src/eko/constants.py
index c01197405..fba786f46 100644
--- a/src/eko/constants.py
+++ b/src/eko/constants.py
@@ -84,6 +84,8 @@ def uplike_flavors(nf):
     nu : int
 
     """
+    if nf > 6:
+        raise NotImplementedError("Selected nf is not implemented")
     nu = nf // 2
     return nu
 

From 16d69f24a33e27ad57e451d29d1d8ecf3d992c0c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Tue, 9 Jan 2024 19:36:44 +0100
Subject: [PATCH 03/22] Implement varying nl in the coupling running

---
 src/eko/constants.py                   | 22 ++++++++++++++++
 src/eko/couplings.py                   | 36 ++++++++++++++++++--------
 src/eko/evolution_operator/__init__.py | 28 ++------------------
 3 files changed, 49 insertions(+), 37 deletions(-)

diff --git a/src/eko/constants.py b/src/eko/constants.py
index fba786f46..1e00d0c62 100644
--- a/src/eko/constants.py
+++ b/src/eko/constants.py
@@ -25,6 +25,9 @@
 Defaults to :math:`\frac{N_C^2-1}{2N_C} = 4/3`.
 """
 
+MTAU = 1.777
+"""Mass of the tau."""
+
 eu2 = 4.0 / 9
 """Up quarks charge squared."""
 
@@ -114,3 +117,22 @@ def charge_combinations(nf):
     vde2m = nd / nf * (eu2 - ed2)
     e2delta = vde2m - vue2m + e2avg
     return e2avg, vue2m, vde2m, e2delta
+
+
+@nb.njit(cache=True)
+def lepton_number(q2):
+    """Compute the number of up flavors.
+
+    at the moment we never go below 1 GeV so we don't need muon and electron
+
+    Parameters
+    ----------
+    q : float
+        energy
+
+    Returns
+    -------
+    nl : int
+
+    """
+    return 3 if q2 > MTAU**2 else 2
diff --git a/src/eko/couplings.py b/src/eko/couplings.py
index b0d39a070..bd6e658db 100644
--- a/src/eko/couplings.py
+++ b/src/eko/couplings.py
@@ -277,7 +277,7 @@ def expanded_qed(ref, order, beta0, b_vec, lmu):
 
 @nb.njit(cache=True)
 def couplings_expanded_alphaem_running(
-    order, couplings_ref, nf, scale_from, scale_to, decoupled_running
+    order, couplings_ref, nf, nl, scale_from, scale_to, decoupled_running
 ):
     r"""Compute coupled expanded expression of the couplings for running alphaem.
 
@@ -308,8 +308,8 @@ def couplings_expanded_alphaem_running(
     beta0_qcd = beta_qcd((2, 0), nf)
     b_vec_qcd = [b_qcd((i + 2, 0), nf) for i in range(order[0])]
     res_as = expanded_qcd(couplings_ref[0], order[0], beta0_qcd, b_vec_qcd, lmu)
-    beta0_qed = beta_qed((0, 2), nf)
-    b_vec_qed = [b_qed((0, i + 2), nf) for i in range(order[1])]
+    beta0_qed = beta_qed((0, 2), nf, nl)
+    b_vec_qed = [b_qed((0, i + 2), nf, nl) for i in range(order[1])]
     res_aem = expanded_qed(couplings_ref[1], order[1], beta0_qed, b_vec_qed, lmu)
     # if order[0] >= 1 and order[1] >= 1:
     # order[0] is always >=1
@@ -521,7 +521,7 @@ def rge(_t, a, b_vec):
         )
         return res.y[0][-1]
 
-    def compute_exact_alphaem_running(self, a_ref, nf, scale_from, scale_to):
+    def compute_exact_alphaem_running(self, a_ref, nf, nl, scale_from, scale_to):
         """Compute couplings via |RGE| with running alphaem.
 
         Parameters
@@ -576,12 +576,12 @@ def compute_exact_alphaem_running(self, a_ref, nf, scale_from, scale_to):
             # while aem is constant
             return np.array([rge_qcd, a_ref[1]])
         if self.order[1] >= 1:
-            beta_qed_vec = [beta_qed((0, 2), nf)]
+            beta_qed_vec = [beta_qed((0, 2), nf, nl)]
             if not self.decoupled_running:
                 beta_qcd_mix = beta_qcd((2, 1), nf)
-                beta_qed_mix = beta_qed((1, 2), nf)  # order[0] is always at least 1
+                beta_qed_mix = beta_qed((1, 2), nf, nl)  # order[0] is always at least 1
             if self.order[1] >= 2:
-                beta_qed_vec.append(beta_qed((0, 3), nf))
+                beta_qed_vec.append(beta_qed((0, 3), nf, nl))
 
         # integration kernel
         def rge(_t, a, beta_qcd_vec, beta_qcd_mix, beta_qed_vec, beta_qed_mix):
@@ -656,7 +656,7 @@ def compute_exact_fixed_alphaem(self, a_ref, nf, scale_from, scale_to):
         )
         return np.array([rge_qcd, a_ref[1]])
 
-    def compute(self, a_ref, nf, scale_from, scale_to):
+    def compute(self, a_ref, nf, nl, scale_from, scale_to):
         """Compute actual couplings.
 
         This is a wrapper in order to pass the computation to the corresponding method
@@ -678,7 +678,7 @@ def compute(self, a_ref, nf, scale_from, scale_to):
         numpy.ndarray
             couplings at target scale :math:`a(Q^2)`
         """
-        key = (float(a_ref[0]), float(a_ref[1]), nf, scale_from, float(scale_to))
+        key = (float(a_ref[0]), float(a_ref[1]), nf, nl, scale_from, float(scale_to))
         try:
             return self.cache[key].copy()
         except KeyError:
@@ -686,7 +686,7 @@ def compute(self, a_ref, nf, scale_from, scale_to):
             if self.method == "exact":
                 if self.alphaem_running:
                     a_new = self.compute_exact_alphaem_running(
-                        a_ref.astype(float), nf, scale_from, scale_to
+                        a_ref.astype(float), nf, nl, scale_from, scale_to
                     )
                 else:
                     a_new = self.compute_exact_fixed_alphaem(
@@ -698,6 +698,7 @@ def compute(self, a_ref, nf, scale_from, scale_to):
                         self.order,
                         a_ref.astype(float),
                         nf,
+                        nl,
                         scale_from,
                         float(scale_to),
                         self.decoupled_running,
@@ -737,7 +738,20 @@ def a(
         for k, seg in enumerate(path):
             # skip a very short segment, but keep the matching
             if not np.isclose(seg.origin, seg.target):
-                new_a = self.compute(final_a, seg.nf, seg.origin, seg.target)
+                nli = constants.lepton_number(seg.origin)
+                nlf = constants.lepton_number(seg.target)
+                if nli != nlf:
+                    # it means that MTAU is between origin and target:
+                    # first we evolve from origin to MTAU with nli leptons
+                    a_tmp = self.compute(
+                        final_a, seg.nf, nli, seg.origin, constants.MTAU**2
+                    )
+                    # then from MTAU to target with nlf leptons
+                    new_a = self.compute(
+                        a_tmp, seg.nf, nlf, constants.MTAU**2, seg.target
+                    )
+                else:
+                    new_a = self.compute(final_a, seg.nf, nli, seg.origin, seg.target)
             else:
                 new_a = final_a
             # apply matching conditions: see hep-ph/9706430
diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index 29e67c19c..8092d31a0 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -704,38 +704,14 @@ def compute_aem_list(self):
             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]
-            aem0 = self.a_em[0]
-            aem1 = self.a_em[1]
-            q2ref = self.managers["couplings"].mu2_ref
-            delta_from = abs(self.q2_from - q2ref)
-            delta_to = abs(self.q2_to - q2ref)
-            # I compute the values in aem_list starting from the mu2
-            # that is closer to mu_ref.
-            if delta_from > delta_to:
-                a_start = np.array([as1, aem1])
-                mu2_start = self.q2_to
-            else:
-                a_start = np.array([as0, aem0])
-                mu2_start = self.q2_from
             couplings = self.managers["couplings"]
             mu2_steps = utils.geomspace(self.q2_from, self.q2_to, 1 + ev_op_iterations)
             mu2_l = mu2_steps[0]
-            as_list = np.array(
-                [
-                    couplings.compute(
-                        a_ref=a_start, nf=self.nf, scale_from=mu2_start, scale_to=mu2
-                    )[0]
-                    for mu2 in mu2_steps
-                ]
-            )
+            as_list = np.array([couplings.a_s(scale_to=mu2)[0] for mu2 in mu2_steps])
             a_half = np.zeros((ev_op_iterations, 2))
             for step, mu2_h in enumerate(mu2_steps[1:]):
                 mu2_half = (mu2_h + mu2_l) / 2.0
-                a_s, aem = couplings.compute(
-                    a_ref=a_start, nf=self.nf, scale_from=mu2_start, scale_to=mu2_half
-                )
+                a_s, aem = couplings.a(scale_to=mu2_half)
                 a_half[step] = [a_s, aem]
                 mu2_l = mu2_h
         return as_list, a_half

From 3bc4fb771024ab0aaf0f9ead64d394d1de9d5eb0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Tue, 9 Jan 2024 20:40:40 +0100
Subject: [PATCH 04/22] Fix tests

---
 src/eko/beta.py        |  8 ++++----
 src/eko/couplings.py   |  2 +-
 tests/eko/test_beta.py | 19 ++++++++++++-------
 3 files changed, 17 insertions(+), 12 deletions(-)

diff --git a/src/eko/beta.py b/src/eko/beta.py
index 6d8b77ec2..d35af42e2 100644
--- a/src/eko/beta.py
+++ b/src/eko/beta.py
@@ -32,7 +32,7 @@ def beta_qcd_as2(nf):
 
 
 @nb.njit(cache=True)
-def beta_qed_aem2(nf, nl=3):
+def beta_qed_aem2(nf, nl):
     r"""Compute the first coefficient of the QED beta function.
 
     Implements Eq. (7) of :cite:`Surguladze:1996hx`.
@@ -82,7 +82,7 @@ def beta_qcd_as3(nf):
 
 
 @nb.njit(cache=True)
-def beta_qed_aem3(nf, nl=3):
+def beta_qed_aem3(nf, nl):
     r"""Compute the second coefficient of the QED beta function.
 
     Implements Eq. (7) of :cite:`Surguladze:1996hx`.
@@ -244,7 +244,7 @@ def beta_qcd(k, nf):
 
 
 @nb.njit(cache=True)
-def beta_qed(k, nf, nl=3):
+def beta_qed(k, nf, nl):
     r"""Compute value of a beta_qed coefficients.
 
     Parameters
@@ -293,7 +293,7 @@ def b_qcd(k, nf):
 
 
 @nb.njit(cache=True)
-def b_qed(k, nf, nl=3):
+def b_qed(k, nf, nl):
     r"""Compute b_qed coefficient.
 
     Parameters
diff --git a/src/eko/couplings.py b/src/eko/couplings.py
index bd6e658db..bdc653832 100644
--- a/src/eko/couplings.py
+++ b/src/eko/couplings.py
@@ -322,7 +322,7 @@ def couplings_expanded_alphaem_running(
             )
             res_aem += (
                 -couplings_ref[1] ** 2
-                * b_qed((1, 2), nf)
+                * b_qed((1, 2), nf, nl)
                 * np.log(1 + beta0_qed * couplings_ref[0] * lmu)
             )
     return np.array([res_as, res_aem])
diff --git a/tests/eko/test_beta.py b/tests/eko/test_beta.py
index cd64c2b85..8c6be7c5c 100644
--- a/tests/eko/test_beta.py
+++ b/tests/eko/test_beta.py
@@ -28,7 +28,7 @@ def test_beta_aem2():
     """Test first beta function coefficient"""
     # from hep-ph/9803211
     np.testing.assert_approx_equal(
-        beta.beta_qed_aem2(5), -4.0 / 3 * (3 + 3 * (2 * 4 / 9 + 3 * 1 / 9))
+        beta.beta_qed_aem2(5, 3), -4.0 / 3 * (3 + 3 * (2 * 4 / 9 + 3 * 1 / 9))
     )
 
 
@@ -43,7 +43,7 @@ def test_beta_aem3():
     """Test second beta function coefficient"""
     # from hep-ph/9803211
     np.testing.assert_approx_equal(
-        beta.beta_qed_aem3(5), -4.0 * (3 + 3 * (2 * 16 / 81 + 3 * 1 / 81))
+        beta.beta_qed_aem3(5, 3), -4.0 * (3 + 3 * (2 * 16 / 81 + 3 * 1 / 81))
     )
 
 
@@ -66,21 +66,26 @@ def test_beta_as5():
 def test_beta():
     """beta-wrapper"""
     nf = 3
+    nl = 3
     np.testing.assert_allclose(beta.beta_qcd((2, 0), nf), beta.beta_qcd_as2(nf))
     np.testing.assert_allclose(beta.beta_qcd((3, 0), nf), beta.beta_qcd_as3(nf))
     np.testing.assert_allclose(beta.beta_qcd((4, 0), nf), beta.beta_qcd_as4(nf))
     np.testing.assert_allclose(beta.beta_qcd((5, 0), nf), beta.beta_qcd_as5(nf))
     np.testing.assert_allclose(beta.beta_qcd((2, 1), nf), beta.beta_qcd_as2aem1(nf))
-    np.testing.assert_allclose(beta.beta_qed((0, 2), nf), beta.beta_qed_aem2(nf))
-    np.testing.assert_allclose(beta.beta_qed((0, 3), nf), beta.beta_qed_aem3(nf))
-    np.testing.assert_allclose(beta.beta_qed((1, 2), nf), beta.beta_qed_aem2as1(nf))
+    np.testing.assert_allclose(
+        beta.beta_qed((0, 2), nf, nl), beta.beta_qed_aem2(nf, nl)
+    )
+    np.testing.assert_allclose(
+        beta.beta_qed((0, 3), nf, nl), beta.beta_qed_aem3(nf, nl)
+    )
+    np.testing.assert_allclose(beta.beta_qed((1, 2), nf, nl), beta.beta_qed_aem2as1(nf))
     with pytest.raises(ValueError):
         beta.beta_qcd((6, 0), 3)
     with pytest.raises(ValueError):
-        beta.beta_qed((0, 4), 3)
+        beta.beta_qed((0, 4), 3, 3)
 
 
 def test_b():
     """b-wrapper"""
     np.testing.assert_allclose(beta.b_qcd((2, 0), 3), 1.0)
-    np.testing.assert_allclose(beta.b_qed((0, 2), 3), 1.0)
+    np.testing.assert_allclose(beta.b_qed((0, 2), 3, 3), 1.0)

From 0f46f5851f99b703a24f8e5f696c1f88f9827dc7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Tue, 9 Jan 2024 21:10:19 +0100
Subject: [PATCH 05/22] Fix test sv

---
 src/eko/evolution_operator/__init__.py           | 8 ++++----
 src/eko/scale_variations/exponentiated.py        | 4 ++--
 tests/eko/scale_variations/test_exponentiated.py | 8 ++++----
 3 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index 8092d31a0..92508f164 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -18,7 +18,7 @@
 import ekore.anomalous_dimensions.unpolarized.time_like as ad_ut
 
 from .. import basis_rotation as br
-from .. import interpolation, mellin
+from .. import constants, interpolation, mellin
 from .. import scale_variations as sv
 from ..kernels import non_singlet as ns
 from ..kernels import non_singlet_qed as qed_ns
@@ -505,7 +505,7 @@ def quad_ker_qed(
         # scale var exponentiated is directly applied on gamma
         if sv_mode == sv.Modes.exponentiated:
             gamma_s = sv.exponentiated.gamma_variation_qed(
-                gamma_s, order, nf, L, alphaem_running
+                gamma_s, order, nf, constants.lepton_number(mu2_to), L, alphaem_running
             )
         ker = qed_s.dispatcher(
             order,
@@ -533,7 +533,7 @@ def quad_ker_qed(
         # scale var exponentiated is directly applied on gamma
         if sv_mode == sv.Modes.exponentiated:
             gamma_v = sv.exponentiated.gamma_variation_qed(
-                gamma_v, order, nf, L, alphaem_running
+                gamma_v, order, nf, constants.lepton_number(mu2_to), L, alphaem_running
             )
         ker = qed_v.dispatcher(
             order,
@@ -558,7 +558,7 @@ def quad_ker_qed(
         # scale var exponentiated is directly applied on gamma
         if sv_mode == sv.Modes.exponentiated:
             gamma_ns = sv.exponentiated.gamma_variation_qed(
-                gamma_ns, order, nf, L, alphaem_running
+                gamma_ns, order, nf, constants.lepton_number(mu2_to), L, alphaem_running
             )
         ker = qed_ns.dispatcher(
             order,
diff --git a/src/eko/scale_variations/exponentiated.py b/src/eko/scale_variations/exponentiated.py
index 13d21c34e..491c71395 100644
--- a/src/eko/scale_variations/exponentiated.py
+++ b/src/eko/scale_variations/exponentiated.py
@@ -46,7 +46,7 @@ def gamma_variation(gamma, order, nf, L):
 
 
 @nb.njit(cache=True)
-def gamma_variation_qed(gamma, order, nf, L, alphaem_running):
+def gamma_variation_qed(gamma, order, nf, nl, L, alphaem_running):
     """Adjust the anomalous dimensions with the scale variations.
 
     Parameters
@@ -71,6 +71,6 @@ def gamma_variation_qed(gamma, order, nf, L, alphaem_running):
     gamma[1:, 0] = gamma_variation(gamma[1:, 0], order, nf, L)
     if alphaem_running:
         if order[1] >= 2:
-            beta0qed = beta.beta_qed((0, 2), nf)
+            beta0qed = beta.beta_qed((0, 2), nf, nl)
             gamma[0, 2] += beta0qed * gamma[0, 1] * L
         return gamma
diff --git a/tests/eko/scale_variations/test_exponentiated.py b/tests/eko/scale_variations/test_exponentiated.py
index 79cc6b386..ee76f9ed9 100644
--- a/tests/eko/scale_variations/test_exponentiated.py
+++ b/tests/eko/scale_variations/test_exponentiated.py
@@ -39,13 +39,13 @@ def test_gamma_ns_qed_fact():
             [0.15, 0.0, 0.0, 0.0, 0.0],
         ]
     ).transpose()
-    gamma_ns_as1aem1_0 = gamma_variation_qed(gamma_ns.copy(), (1, 1), 3, 0, True)
+    gamma_ns_as1aem1_0 = gamma_variation_qed(gamma_ns.copy(), (1, 1), 3, 3, 0, True)
     np.testing.assert_allclose(gamma_ns_as1aem1_0, gamma_ns)
-    gamma_ns_as1aem1_1 = gamma_variation_qed(gamma_ns.copy(), (1, 1), 3, -1, True)
+    gamma_ns_as1aem1_1 = gamma_variation_qed(gamma_ns.copy(), (1, 1), 3, 3, -1, True)
     np.testing.assert_allclose(gamma_ns_as1aem1_1, gamma_ns)
-    gamma_ns_as2aem2_1 = gamma_variation_qed(gamma_ns.copy(), (2, 2), 3, -1, True)
+    gamma_ns_as2aem2_1 = gamma_variation_qed(gamma_ns.copy(), (2, 2), 3, 3, -1, True)
     assert gamma_ns_as2aem2_1[2, 0] < gamma_ns[2, 0]
     assert gamma_ns_as2aem2_1[0, 2] > gamma_ns[0, 2]  # beta0qed < 0
-    gamma_ns_as4aem2_1 = gamma_variation_qed(gamma_ns.copy(), (4, 2), 3, -1, True)
+    gamma_ns_as4aem2_1 = gamma_variation_qed(gamma_ns.copy(), (4, 2), 3, 3, -1, True)
     gamma_ns_N3LO_1 = gamma_variation(gamma_ns.copy()[1:, 0], (4, 0), 3, -1)
     np.testing.assert_allclose(gamma_ns_as4aem2_1[1:, 0], gamma_ns_N3LO_1)

From 3628b3fdb90f336632818cf09d60a5dd0e4780ef Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Tue, 9 Jan 2024 21:20:11 +0100
Subject: [PATCH 06/22] Fix typo in evolution_operator/init

---
 src/eko/evolution_operator/__init__.py    | 2 +-
 tests/eko/evolution_operator/test_init.py | 6 ++++++
 2 files changed, 7 insertions(+), 1 deletion(-)

diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index 92508f164..aa67ecbbe 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -707,7 +707,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]
-            as_list = np.array([couplings.a_s(scale_to=mu2)[0] for mu2 in mu2_steps])
+            as_list = np.array([couplings.a_s(scale_to=mu2) for mu2 in mu2_steps])
             a_half = np.zeros((ev_op_iterations, 2))
             for step, mu2_h in enumerate(mu2_steps[1:]):
                 mu2_half = (mu2_h + mu2_l) / 2.0
diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py
index 47d5e7001..e1a3ac705 100644
--- a/tests/eko/evolution_operator/test_init.py
+++ b/tests/eko/evolution_operator/test_init.py
@@ -197,6 +197,12 @@ def __init__(self):
     def a(self, scale_to=None, nf_to=None):
         return (0.1, 0.01)
 
+    def a_s(self, scale_to=None, nf_to=None):
+        return self.a(scale_to, nf_to)[0]
+
+    def a_em(self, scale_to=None, nf_to=None):
+        return self.a(scale_to, nf_to)[1]
+
     def compute(self, a_ref, nf, scale_from, scale_to):
         return a_ref
 

From bc956e25c3c7cab554088f21ccf20c6ddd9cc7c1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Tue, 9 Jan 2024 22:05:24 +0100
Subject: [PATCH 07/22] Increase rtol to make bench pass

---
 benchmarks/eko/benchmark_msbar_evolution.py | 2 +-
 benchmarks/eko/benchmark_strong_coupling.py | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/benchmarks/eko/benchmark_msbar_evolution.py b/benchmarks/eko/benchmark_msbar_evolution.py
index 6286e15ad..1e7349ff7 100644
--- a/benchmarks/eko/benchmark_msbar_evolution.py
+++ b/benchmarks/eko/benchmark_msbar_evolution.py
@@ -174,7 +174,7 @@ def benchmark_APFEL_msbar_evolution(
                     )
                 # check myself to APFEL
                 np.testing.assert_allclose(
-                    apfel_vals, np.sqrt(np.array(my_vals)), rtol=2.3e-3
+                    apfel_vals, np.sqrt(np.array(my_vals)), rtol=2.6e-3
                 )
 
     def benchmark_APFEL_msbar_solution(
diff --git a/benchmarks/eko/benchmark_strong_coupling.py b/benchmarks/eko/benchmark_strong_coupling.py
index d13a264bd..911d529c2 100644
--- a/benchmarks/eko/benchmark_strong_coupling.py
+++ b/benchmarks/eko/benchmark_strong_coupling.py
@@ -322,7 +322,7 @@ def benchmark_APFEL_vfns(self):
                 # print(apfel_vals_cur)
                 np.testing.assert_allclose(apfel_vals, np.array(apfel_vals_cur))
             # check myself to APFEL
-            np.testing.assert_allclose(apfel_vals, np.array(my_vals))
+            np.testing.assert_allclose(apfel_vals, np.array(my_vals), rtol=3e-3)
 
     def benchmark_APFEL_vfns_fact_to_ren(self):
         Q2s = [

From 88c4dcc2f816e26a653496a7b283c50789c7c6ca Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?=
 <57321974+niclaurenti@users.noreply.github.com>
Date: Wed, 10 Jan 2024 12:40:17 +0100
Subject: [PATCH 08/22] Update src/eko/constants.py

Co-authored-by: Felix Hekhorn <felixhekhorn@users.noreply.github.com>
---
 src/eko/constants.py | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/eko/constants.py b/src/eko/constants.py
index 1e00d0c62..b208e4182 100644
--- a/src/eko/constants.py
+++ b/src/eko/constants.py
@@ -121,18 +121,19 @@ def charge_combinations(nf):
 
 @nb.njit(cache=True)
 def lepton_number(q2):
-    """Compute the number of up flavors.
+    """Compute the number of leptons.
 
-    at the moment we never go below 1 GeV so we don't need muon and electron
+    Note: at the moment we never go below 1 GeV so we don't need muon and electron.
 
     Parameters
     ----------
-    q : float
-        energy
+    q2 : float
+        scale
 
     Returns
     -------
-    nl : int
+    int :
+       Number of leptons
 
     """
     return 3 if q2 > MTAU**2 else 2

From 6bcacc4698ad7811de006fb4001f24c80b95e32c Mon Sep 17 00:00:00 2001
From: niccolo <niclaurenti@gmail.com>
Date: Wed, 10 Jan 2024 14:21:09 +0100
Subject: [PATCH 09/22] Update docstring

---
 src/eko/couplings.py                      | 6 ++++++
 src/eko/scale_variations/exponentiated.py | 2 ++
 2 files changed, 8 insertions(+)

diff --git a/src/eko/couplings.py b/src/eko/couplings.py
index bdc653832..8be9f770d 100644
--- a/src/eko/couplings.py
+++ b/src/eko/couplings.py
@@ -291,6 +291,8 @@ def couplings_expanded_alphaem_running(
         reference alpha_s and alpha
     nf : int
         value of nf for computing the couplings
+    nl : int
+        number of leptons partecipating to alphaem running
     scale_from : float
         reference scale
     scale_to : float
@@ -530,6 +532,8 @@ def compute_exact_alphaem_running(self, a_ref, nf, nl, scale_from, scale_to):
             reference alpha_s and alpha
         nf : int
             value of nf for computing alpha_i
+        nl : int
+            number of leptons partecipating to alphaem running
         scale_from : float
             reference scale
         scale_to : float
@@ -668,6 +672,8 @@ def compute(self, a_ref, nf, nl, scale_from, scale_to):
             reference a
         nf : int
             value of nf for computing alpha
+        nl : int
+            number of leptons partecipating to alphaem running
         scale_from : float
             reference scale
         scale_to : float
diff --git a/src/eko/scale_variations/exponentiated.py b/src/eko/scale_variations/exponentiated.py
index 491c71395..7efb88e79 100644
--- a/src/eko/scale_variations/exponentiated.py
+++ b/src/eko/scale_variations/exponentiated.py
@@ -57,6 +57,8 @@ def gamma_variation_qed(gamma, order, nf, nl, L, alphaem_running):
         perturbation order
     nf : int
         number of active flavors
+    nl : int
+        number of leptons partecipating to alphaem running
     L : float
         logarithmic ratio of factorization and renormalization scale
 

From ce868a10ebfed23c2b63f8cf003766c1429f6f07 Mon Sep 17 00:00:00 2001
From: niccolo <niclaurenti@gmail.com>
Date: Wed, 10 Jan 2024 14:47:11 +0100
Subject: [PATCH 10/22] Move lepton_number in matchings

---
 src/eko/constants.py                   | 19 -------------------
 src/eko/evolution_operator/__init__.py |  9 +++++----
 src/eko/matchings.py                   | 18 ++++++++++++++++++
 3 files changed, 23 insertions(+), 23 deletions(-)

diff --git a/src/eko/constants.py b/src/eko/constants.py
index 1e00d0c62..55dff40a2 100644
--- a/src/eko/constants.py
+++ b/src/eko/constants.py
@@ -117,22 +117,3 @@ def charge_combinations(nf):
     vde2m = nd / nf * (eu2 - ed2)
     e2delta = vde2m - vue2m + e2avg
     return e2avg, vue2m, vde2m, e2delta
-
-
-@nb.njit(cache=True)
-def lepton_number(q2):
-    """Compute the number of up flavors.
-
-    at the moment we never go below 1 GeV so we don't need muon and electron
-
-    Parameters
-    ----------
-    q : float
-        energy
-
-    Returns
-    -------
-    nl : int
-
-    """
-    return 3 if q2 > MTAU**2 else 2
diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index aa67ecbbe..12132a091 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -18,7 +18,8 @@
 import ekore.anomalous_dimensions.unpolarized.time_like as ad_ut
 
 from .. import basis_rotation as br
-from .. import constants, interpolation, mellin
+from .. import interpolation, mellin
+from ..matchings import lepton_number
 from .. import scale_variations as sv
 from ..kernels import non_singlet as ns
 from ..kernels import non_singlet_qed as qed_ns
@@ -505,7 +506,7 @@ def quad_ker_qed(
         # scale var exponentiated is directly applied on gamma
         if sv_mode == sv.Modes.exponentiated:
             gamma_s = sv.exponentiated.gamma_variation_qed(
-                gamma_s, order, nf, constants.lepton_number(mu2_to), L, alphaem_running
+                gamma_s, order, nf, lepton_number(mu2_to), L, alphaem_running
             )
         ker = qed_s.dispatcher(
             order,
@@ -533,7 +534,7 @@ def quad_ker_qed(
         # scale var exponentiated is directly applied on gamma
         if sv_mode == sv.Modes.exponentiated:
             gamma_v = sv.exponentiated.gamma_variation_qed(
-                gamma_v, order, nf, constants.lepton_number(mu2_to), L, alphaem_running
+                gamma_v, order, nf, lepton_number(mu2_to), L, alphaem_running
             )
         ker = qed_v.dispatcher(
             order,
@@ -558,7 +559,7 @@ def quad_ker_qed(
         # scale var exponentiated is directly applied on gamma
         if sv_mode == sv.Modes.exponentiated:
             gamma_ns = sv.exponentiated.gamma_variation_qed(
-                gamma_ns, order, nf, constants.lepton_number(mu2_to), L, alphaem_running
+                gamma_ns, order, nf, lepton_number(mu2_to), L, alphaem_running
             )
         ker = qed_ns.dispatcher(
             order,
diff --git a/src/eko/matchings.py b/src/eko/matchings.py
index 0dfcae24f..721727b9f 100644
--- a/src/eko/matchings.py
+++ b/src/eko/matchings.py
@@ -8,6 +8,7 @@
 from .io.types import EvolutionPoint as EPoint
 from .io.types import FlavorIndex, FlavorsNumber, SquaredScale
 from .quantities.heavy_quarks import MatchingScales
+from .constants import MTAU
 
 logger = logging.getLogger(__name__)
 
@@ -209,3 +210,20 @@ def is_downward_path(path: Path) -> bool:
 def flavor_shift(is_downward: bool) -> int:
     """Determine the shift to number of light flavors."""
     return 4 if is_downward else 3
+
+def lepton_number(q2):
+    """Compute the number of up flavors.
+
+    at the moment we never go below 1 GeV so we don't need muon and electron
+
+    Parameters
+    ----------
+    q : float
+        energy
+
+    Returns
+    -------
+    nl : int
+
+    """
+    return 3 if q2 > MTAU**2 else 2

From c493ecc6aad895332065b86938bb0a5625027b0d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Wed, 10 Jan 2024 15:19:28 +0100
Subject: [PATCH 11/22] Add check

---
 benchmarks/eko/benchmark_msbar_evolution.py | 2 +-
 benchmarks/eko/benchmark_strong_coupling.py | 2 +-
 src/eko/couplings.py                        | 2 +-
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/benchmarks/eko/benchmark_msbar_evolution.py b/benchmarks/eko/benchmark_msbar_evolution.py
index 1e7349ff7..6286e15ad 100644
--- a/benchmarks/eko/benchmark_msbar_evolution.py
+++ b/benchmarks/eko/benchmark_msbar_evolution.py
@@ -174,7 +174,7 @@ def benchmark_APFEL_msbar_evolution(
                     )
                 # check myself to APFEL
                 np.testing.assert_allclose(
-                    apfel_vals, np.sqrt(np.array(my_vals)), rtol=2.6e-3
+                    apfel_vals, np.sqrt(np.array(my_vals)), rtol=2.3e-3
                 )
 
     def benchmark_APFEL_msbar_solution(
diff --git a/benchmarks/eko/benchmark_strong_coupling.py b/benchmarks/eko/benchmark_strong_coupling.py
index 911d529c2..d13a264bd 100644
--- a/benchmarks/eko/benchmark_strong_coupling.py
+++ b/benchmarks/eko/benchmark_strong_coupling.py
@@ -322,7 +322,7 @@ def benchmark_APFEL_vfns(self):
                 # print(apfel_vals_cur)
                 np.testing.assert_allclose(apfel_vals, np.array(apfel_vals_cur))
             # check myself to APFEL
-            np.testing.assert_allclose(apfel_vals, np.array(my_vals), rtol=3e-3)
+            np.testing.assert_allclose(apfel_vals, np.array(my_vals))
 
     def benchmark_APFEL_vfns_fact_to_ren(self):
         Q2s = [
diff --git a/src/eko/couplings.py b/src/eko/couplings.py
index bdc653832..b437e88ea 100644
--- a/src/eko/couplings.py
+++ b/src/eko/couplings.py
@@ -740,7 +740,7 @@ def a(
             if not np.isclose(seg.origin, seg.target):
                 nli = constants.lepton_number(seg.origin)
                 nlf = constants.lepton_number(seg.target)
-                if nli != nlf:
+                if self.order[1] != 0 and nli != nlf:
                     # it means that MTAU is between origin and target:
                     # first we evolve from origin to MTAU with nli leptons
                     a_tmp = self.compute(

From 1999d92e71a8ba687f1f611af237fcd88f96e5b4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Wed, 10 Jan 2024 15:54:41 +0100
Subject: [PATCH 12/22] Call isort on src/eko/evolution_operator/__init__.py

---
 src/eko/evolution_operator/__init__.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index 12132a091..1e132c50d 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -19,7 +19,6 @@
 
 from .. import basis_rotation as br
 from .. import interpolation, mellin
-from ..matchings import lepton_number
 from .. import scale_variations as sv
 from ..kernels import non_singlet as ns
 from ..kernels import non_singlet_qed as qed_ns
@@ -27,7 +26,7 @@
 from ..kernels import singlet_qed as qed_s
 from ..kernels import utils
 from ..kernels import valence_qed as qed_v
-from ..matchings import Segment
+from ..matchings import Segment, lepton_number
 from ..member import OpMember
 
 logger = logging.getLogger(__name__)

From c4951c73eefdb3372a3c20975205e16fa5b86317 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?= <niclaurenti@gmail.com>
Date: Wed, 10 Jan 2024 16:17:01 +0100
Subject: [PATCH 13/22] Add numba decorator

---
 src/eko/matchings.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/eko/matchings.py b/src/eko/matchings.py
index 74a2aa478..4e529bb2a 100644
--- a/src/eko/matchings.py
+++ b/src/eko/matchings.py
@@ -3,6 +3,7 @@
 from dataclasses import dataclass
 from typing import List, Union
 
+import numba as nb
 import numpy as np
 
 from .constants import MTAU
@@ -212,6 +213,7 @@ def flavor_shift(is_downward: bool) -> int:
     return 4 if is_downward else 3
 
 
+@nb.njit(cache=True)
 def lepton_number(q2):
     """Compute the number of leptons.
 

From e561cfde2fa22144387fa3f23c5aaa2efb8357ef Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?=
 <57321974+niclaurenti@users.noreply.github.com>
Date: Wed, 10 Jan 2024 16:30:43 +0100
Subject: [PATCH 14/22] Update src/eko/matchings.py

Co-authored-by: Giacomo Magni <39065935+giacomomagni@users.noreply.github.com>
---
 src/eko/matchings.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/eko/matchings.py b/src/eko/matchings.py
index 4e529bb2a..9ef6b8f9e 100644
--- a/src/eko/matchings.py
+++ b/src/eko/matchings.py
@@ -217,7 +217,7 @@ def flavor_shift(is_downward: bool) -> int:
 def lepton_number(q2):
     """Compute the number of leptons.
 
-    Note: at the moment we never go below 1 GeV so we don't need muon and electron.
+    Note: muons and electrons are always massless as for up, down and strange.
 
     Parameters
     ----------

From 75bfca0dda478488ee65ee8f2e4f6db027594955 Mon Sep 17 00:00:00 2001
From: niccolo <niclaurenti@gmail.com>
Date: Wed, 10 Jan 2024 16:50:06 +0100
Subject: [PATCH 15/22] Add docstrings

---
 src/eko/beta.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/src/eko/beta.py b/src/eko/beta.py
index d35af42e2..5f8f88193 100644
--- a/src/eko/beta.py
+++ b/src/eko/beta.py
@@ -41,6 +41,8 @@ def beta_qed_aem2(nf, nl):
     ----------
         nf : int
             number of active flavors
+        nl : int
+            number of leptons partecipating to alphaem running
 
     Returns
     -------
@@ -91,6 +93,8 @@ def beta_qed_aem3(nf, nl):
     ----------
         nf : int
             number of active flavors
+        nl : int
+            number of leptons partecipating to alphaem running
 
     Returns
     -------
@@ -253,6 +257,8 @@ def beta_qed(k, nf, nl):
             perturbative order
         nf : int
             number of active flavors
+        nl : int
+            number of leptons partecipating to alphaem running
 
     Returns
     -------
@@ -302,6 +308,8 @@ def b_qed(k, nf, nl):
             perturbative order
         nf : int
             number of active flavors
+        nl : int
+            number of leptons partecipating to alphaem running
 
     Returns
     -------

From b623f974fc180c3f2e8efb9d74c17d16b8798d84 Mon Sep 17 00:00:00 2001
From: niccolo <niclaurenti@gmail.com>
Date: Wed, 10 Jan 2024 18:18:05 +0100
Subject: [PATCH 16/22] Force argument to self.nf

---
 src/eko/evolution_operator/__init__.py | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index 1e132c50d..2142345ef 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -707,13 +707,17 @@ 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]
-            as_list = np.array([couplings.a_s(scale_to=mu2) for mu2 in mu2_steps])
+            as_list = np.array(
+                [couplings.a_s(scale_to=mu2, nf_to=self.nf) for mu2 in mu2_steps]
+            )
             a_half = np.zeros((ev_op_iterations, 2))
             for step, mu2_h in enumerate(mu2_steps[1:]):
                 mu2_half = (mu2_h + mu2_l) / 2.0
-                a_s, aem = couplings.a(scale_to=mu2_half)
+                a_s, aem = couplings.a(scale_to=mu2_half, nf_to=self.nf)
                 a_half[step] = [a_s, aem]
                 mu2_l = mu2_h
+        print(as_list.tolist())
+        print(a_half.tolist())
         return as_list, a_half
 
     @property

From 0e854446ffcbe600ad9ec08b52d9538b4e83a60b Mon Sep 17 00:00:00 2001
From: niccolo <niclaurenti@gmail.com>
Date: Thu, 11 Jan 2024 10:05:22 +0100
Subject: [PATCH 17/22] Remove debugging print

---
 src/eko/evolution_operator/__init__.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py
index 2142345ef..07a460b83 100644
--- a/src/eko/evolution_operator/__init__.py
+++ b/src/eko/evolution_operator/__init__.py
@@ -716,8 +716,6 @@ def compute_aem_list(self):
                 a_s, aem = couplings.a(scale_to=mu2_half, nf_to=self.nf)
                 a_half[step] = [a_s, aem]
                 mu2_l = mu2_h
-        print(as_list.tolist())
-        print(a_half.tolist())
         return as_list, a_half
 
     @property

From efd45623372608304d2eeebfbbe94fe030d7a1dd Mon Sep 17 00:00:00 2001
From: niccolo <niclaurenti@gmail.com>
Date: Thu, 11 Jan 2024 11:54:24 +0100
Subject: [PATCH 18/22] Add bench for alphaem

---
 benchmarks/eko/benchmark_alphaem.py | 144 ++++++++++++++++++++++++++++
 1 file changed, 144 insertions(+)
 create mode 100644 benchmarks/eko/benchmark_alphaem.py

diff --git a/benchmarks/eko/benchmark_alphaem.py b/benchmarks/eko/benchmark_alphaem.py
new file mode 100644
index 000000000..adfbbe581
--- /dev/null
+++ b/benchmarks/eko/benchmark_alphaem.py
@@ -0,0 +1,144 @@
+"""This module benchmarks alpha_em against alphaQED23."""
+
+import numpy as np
+import pytest
+
+from eko import matchings
+from eko.beta import beta_qcd
+from eko.couplings import Couplings
+from eko.io.runcards import TheoryCard
+from eko.quantities.couplings import CouplingEvolutionMethod, CouplingsInfo
+from eko.quantities.heavy_quarks import QuarkMassScheme
+
+
+@pytest.mark.isolated
+class BenchmarkCouplings:
+    def test_alphaQED_high(self):
+        """testing aginst alphaQED23 for high Q values"""
+        alphaQED23 = np.array(
+            [
+                0.0075683,
+                0.0075867,
+                0.0076054,
+                0.0076245,
+                0.0076437,
+                0.0076631,
+                0.0076827,
+                0.0077024,
+                0.0077222,
+                0.0077421,
+                0.0077621,
+            ]
+        )
+        qvec = np.geomspace(10, 100, 11)
+        couplings = CouplingsInfo.from_dict(
+            dict(
+                alphas=0.118,
+                alphaem=7.7553e-03,
+                scale=91.2,
+                num_flavs_ref=5,
+                max_num_flavs=5,
+                em_running=True,
+            )
+        )
+        eko_alpha = Couplings(
+            couplings,
+            (3, 2),
+            method=CouplingEvolutionMethod.EXACT,
+            masses=[m**2 for m in [1.51, 4.92, 172.5]],
+            hqm_scheme=QuarkMassScheme.POLE,
+            thresholds_ratios=[1.0, 1.0, np.inf],
+        )
+        alpha_eko = np.array([eko_alpha.a_em(q**2) * 4 * np.pi for q in qvec])
+
+        np.testing.assert_allclose(alphaQED23, alpha_eko, rtol=1.8e-4)
+
+    def test_alphaQED_low(self):
+        """testing aginst alphaQED23 for low Q values: they are close but not identical"""
+        alphaQED23 = np.array(
+            [
+                0.0074192,
+                0.007431,
+                0.0074434,
+                0.0074565,
+                0.0074702,
+                0.0074847,
+                0.0075,
+                0.0075161,
+                0.0075329,
+                0.0075503,
+                0.0075683,
+            ]
+        )
+        qvec = np.geomspace(1, 10, 11)
+        couplings = CouplingsInfo.from_dict(
+            dict(
+                alphas=0.118,
+                alphaem=7.7553e-03,
+                scale=91.2,
+                num_flavs_ref=5,
+                max_num_flavs=5,
+                em_running=True,
+            )
+        )
+        eko_alpha = Couplings(
+            couplings,
+            (3, 2),
+            method=CouplingEvolutionMethod.EXACT,
+            masses=[m**2 for m in [1.51, 4.92, 172.5]],
+            hqm_scheme=QuarkMassScheme.POLE,
+            thresholds_ratios=[1.0, 1.0, np.inf],
+        )
+        alpha_eko = np.array([eko_alpha.a_em(q**2) * 4 * np.pi for q in qvec])
+
+        np.testing.assert_allclose(alphaQED23, alpha_eko, rtol=3.2e-3)
+
+    def test_validphys(self):
+        """testing aginst validphys"""
+        alpha_vp = np.array(
+            [
+                0.007408255774054356,
+                0.007425240094018394,
+                0.007449051094996458,
+                0.007476301027742958,
+                0.007503751810862984,
+                0.007532299008699658,
+                0.007561621958607614,
+                0.007591174885612722,
+                0.007620960508577136,
+                0.00765098158940323,
+                0.007681240933888789,
+                0.007711741392602655,
+                0.007742485861781425,
+                0.007773477284247778,
+                0.007804718650351058,
+                0.007836212998930739,
+                0.00786796341830342,
+                0.007899973047274033,
+                0.007932245076171957,
+                0.00796478274791273,
+            ]
+        )
+        qvec = np.geomspace(1, 1000, 20)
+        couplings = CouplingsInfo.from_dict(
+            dict(
+                alphas=0.118,
+                alphaem=7.7553e-03,
+                scale=91.2,
+                num_flavs_ref=5,
+                max_num_flavs=5,
+                em_running=True,
+            )
+        )
+        eko_alpha = Couplings(
+            couplings,
+            (3, 2),
+            method=CouplingEvolutionMethod.EXACT,
+            masses=[m**2 for m in [1.51, 4.92, 172.5]],
+            hqm_scheme=QuarkMassScheme.POLE,
+            thresholds_ratios=[1.0, 1.0, np.inf],
+        )
+        eko_alpha.decoupled_running = True
+        alpha_eko = np.array([eko_alpha.a_em(q**2) * 4 * np.pi for q in qvec])
+
+        np.testing.assert_allclose(alpha_vp, alpha_eko, rtol=5e-6)

From e52d0cd73056c316a7b262f230b0fabc3bea3cb1 Mon Sep 17 00:00:00 2001
From: niccolo <niclaurenti@gmail.com>
Date: Thu, 11 Jan 2024 12:55:43 +0100
Subject: [PATCH 19/22] Add cases in test_beta

---
 tests/eko/test_beta.py | 52 ++++++++++++++++++++++++++----------------
 1 file changed, 32 insertions(+), 20 deletions(-)

diff --git a/tests/eko/test_beta.py b/tests/eko/test_beta.py
index 8c6be7c5c..5f1dbd53d 100644
--- a/tests/eko/test_beta.py
+++ b/tests/eko/test_beta.py
@@ -6,7 +6,7 @@
 import pytest
 
 from eko import beta
-from eko.constants import zeta3
+from eko.constants import uplike_flavors, zeta3
 
 
 def _flav_test(function):
@@ -27,9 +27,14 @@ def test_beta_as2():
 def test_beta_aem2():
     """Test first beta function coefficient"""
     # from hep-ph/9803211
-    np.testing.assert_approx_equal(
-        beta.beta_qed_aem2(5, 3), -4.0 / 3 * (3 + 3 * (2 * 4 / 9 + 3 * 1 / 9))
-    )
+    for nf in range(3, 6 + 1):
+        for nl in range(2, 3 + 1):
+            nu = uplike_flavors(nf)
+            nd = nf - nu
+            np.testing.assert_approx_equal(
+                beta.beta_qed_aem2(nf, nl),
+                -4.0 / 3 * (nl + 3 * (nu * 4 / 9 + nd * 1 / 9)),
+            )
 
 
 def test_beta_as3():
@@ -42,9 +47,14 @@ def test_beta_as3():
 def test_beta_aem3():
     """Test second beta function coefficient"""
     # from hep-ph/9803211
-    np.testing.assert_approx_equal(
-        beta.beta_qed_aem3(5, 3), -4.0 * (3 + 3 * (2 * 16 / 81 + 3 * 1 / 81))
-    )
+    for nf in range(3, 6 + 1):
+        for nl in range(2, 3 + 1):
+            nu = uplike_flavors(nf)
+            nd = nf - nu
+            np.testing.assert_approx_equal(
+                beta.beta_qed_aem3(nf, nl),
+                -4.0 * (nl + 3 * (nu * 16 / 81 + nd * 1 / 81)),
+            )
 
 
 def test_beta_as4():
@@ -65,19 +75,21 @@ def test_beta_as5():
 
 def test_beta():
     """beta-wrapper"""
-    nf = 3
-    nl = 3
-    np.testing.assert_allclose(beta.beta_qcd((2, 0), nf), beta.beta_qcd_as2(nf))
-    np.testing.assert_allclose(beta.beta_qcd((3, 0), nf), beta.beta_qcd_as3(nf))
-    np.testing.assert_allclose(beta.beta_qcd((4, 0), nf), beta.beta_qcd_as4(nf))
-    np.testing.assert_allclose(beta.beta_qcd((5, 0), nf), beta.beta_qcd_as5(nf))
-    np.testing.assert_allclose(beta.beta_qcd((2, 1), nf), beta.beta_qcd_as2aem1(nf))
-    np.testing.assert_allclose(
-        beta.beta_qed((0, 2), nf, nl), beta.beta_qed_aem2(nf, nl)
-    )
-    np.testing.assert_allclose(
-        beta.beta_qed((0, 3), nf, nl), beta.beta_qed_aem3(nf, nl)
-    )
+    for nf in range(3, 6 + 1):
+        for nl in range(2, 3 + 1):
+            np.testing.assert_allclose(beta.beta_qcd((2, 0), nf), beta.beta_qcd_as2(nf))
+            np.testing.assert_allclose(beta.beta_qcd((3, 0), nf), beta.beta_qcd_as3(nf))
+            np.testing.assert_allclose(beta.beta_qcd((4, 0), nf), beta.beta_qcd_as4(nf))
+            np.testing.assert_allclose(beta.beta_qcd((5, 0), nf), beta.beta_qcd_as5(nf))
+            np.testing.assert_allclose(
+                beta.beta_qcd((2, 1), nf), beta.beta_qcd_as2aem1(nf)
+            )
+            np.testing.assert_allclose(
+                beta.beta_qed((0, 2), nf, nl), beta.beta_qed_aem2(nf, nl)
+            )
+            np.testing.assert_allclose(
+                beta.beta_qed((0, 3), nf, nl), beta.beta_qed_aem3(nf, nl)
+            )
     np.testing.assert_allclose(beta.beta_qed((1, 2), nf, nl), beta.beta_qed_aem2as1(nf))
     with pytest.raises(ValueError):
         beta.beta_qcd((6, 0), 3)

From 8910ab57622e82b0db5084044705e4eff8d50757 Mon Sep 17 00:00:00 2001
From: niccolo <niclaurenti@gmail.com>
Date: Thu, 11 Jan 2024 14:49:45 +0100
Subject: [PATCH 20/22] add tests

---
 tests/eko/test_beta.py | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

diff --git a/tests/eko/test_beta.py b/tests/eko/test_beta.py
index 5f1dbd53d..682c79f57 100644
--- a/tests/eko/test_beta.py
+++ b/tests/eko/test_beta.py
@@ -101,3 +101,18 @@ def test_b():
     """b-wrapper"""
     np.testing.assert_allclose(beta.b_qcd((2, 0), 3), 1.0)
     np.testing.assert_allclose(beta.b_qed((0, 2), 3, 3), 1.0)
+
+
+def test_zero():
+    """test that beta_qed is zero for nf=nl=0"""
+    for ord in range(2, 3 + 1):
+        np.testing.assert_allclose(beta.beta_qed((0, ord), 0, 0), 0)
+
+
+def test_linearity():
+    """test linearity of QED beta function when nf=0"""
+    for ord in range(2, 3 + 1):
+        for nl in [0, 2]:
+            np.testing.assert_allclose(
+                beta.beta_qed((0, ord), 0, nl) / 2, beta.beta_qed((0, ord), 0, nl / 2)
+            )

From 981af57b49dff20aa5a96f3ea1f99963ea957bf3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Niccol=C3=B2=20Laurenti?=
 <57321974+niclaurenti@users.noreply.github.com>
Date: Thu, 11 Jan 2024 15:05:16 +0100
Subject: [PATCH 21/22] Update benchmarks/eko/benchmark_alphaem.py

Co-authored-by: Felix Hekhorn <felixhekhorn@users.noreply.github.com>
---
 benchmarks/eko/benchmark_alphaem.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/benchmarks/eko/benchmark_alphaem.py b/benchmarks/eko/benchmark_alphaem.py
index adfbbe581..c4d67a1bc 100644
--- a/benchmarks/eko/benchmark_alphaem.py
+++ b/benchmarks/eko/benchmark_alphaem.py
@@ -1,4 +1,7 @@
-"""This module benchmarks alpha_em against alphaQED23."""
+"""This module benchmarks alpha_em against alphaQED23.
+
+alphaQED23 can be obtained from http://www-com.physik.hu-berlin.de/~fjeger/software.html .
+"""
 
 import numpy as np
 import pytest

From 9e23e47e87e642dc15dbef44c00be840913df6dc Mon Sep 17 00:00:00 2001
From: niccolo <niclaurenti@gmail.com>
Date: Fri, 12 Jan 2024 11:19:11 +0100
Subject: [PATCH 22/22] Fix docs of bench_alphaem

---
 benchmarks/eko/benchmark_alphaem.py | 5 +----
 1 file changed, 1 insertion(+), 4 deletions(-)

diff --git a/benchmarks/eko/benchmark_alphaem.py b/benchmarks/eko/benchmark_alphaem.py
index c4d67a1bc..16e56ac4a 100644
--- a/benchmarks/eko/benchmark_alphaem.py
+++ b/benchmarks/eko/benchmark_alphaem.py
@@ -1,4 +1,4 @@
-"""This module benchmarks alpha_em against alphaQED23.
+"""This module benchmarks alpha_em against alphaQED23 and validphys.
 
 alphaQED23 can be obtained from http://www-com.physik.hu-berlin.de/~fjeger/software.html .
 """
@@ -6,10 +6,7 @@
 import numpy as np
 import pytest
 
-from eko import matchings
-from eko.beta import beta_qcd
 from eko.couplings import Couplings
-from eko.io.runcards import TheoryCard
 from eko.quantities.couplings import CouplingEvolutionMethod, CouplingsInfo
 from eko.quantities.heavy_quarks import QuarkMassScheme