Pass nl as an argument in qed beta functions

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):
     r"""Compute the first coefficient of the QED beta function.
 
     Implements Eq. (7) of :cite:`Surguladze:1996hx`.
@@ -41,6 +41,8 @@ def beta_qed_aem2(nf):
     ----------
         nf : int
             number of active flavors
+        nl : int
+            number of leptons partecipating to alphaem running
 
     Returns
     -------
@@ -50,7 +52,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 +84,7 @@ def beta_qcd_as3(nf):
 
 
 @nb.njit(cache=True)
-def beta_qed_aem3(nf):
+def beta_qed_aem3(nf, nl):
     r"""Compute the second coefficient of the QED beta function.
 
     Implements Eq. (7) of :cite:`Surguladze:1996hx`.
@@ -92,6 +93,8 @@ def beta_qed_aem3(nf):
     ----------
         nf : int
             number of active flavors
+        nl : int
+            number of leptons partecipating to alphaem running
 
     Returns
     -------
@@ -101,7 +104,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 +248,7 @@ def beta_qcd(k, nf):
 
 
 @nb.njit(cache=True)
-def beta_qed(k, nf):
+def beta_qed(k, nf, nl):
     r"""Compute value of a beta_qed coefficients.
 
     Parameters
@@ -255,6 +257,8 @@ def beta_qed(k, nf):
             perturbative order
         nf : int
             number of active flavors
+        nl : int
+            number of leptons partecipating to alphaem running
 
     Returns
     -------
@@ -264,9 +268,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 +299,7 @@ def b_qcd(k, nf):
 
 
 @nb.njit(cache=True)
-def b_qed(k, nf):
+def b_qed(k, nf, nl):
     r"""Compute b_qed coefficient.
 
     Parameters
@@ -304,11 +308,13 @@ def b_qed(k, nf):
             perturbative order
         nf : int
             number of active flavors
+        nl : int
+            number of leptons partecipating to alphaem running
 
     Returns
     -------
         b_qed : float
             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)

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."""
 
@@ -84,7 +87,7 @@ def uplike_flavors(nf):
     nu : int
 
     """
-    if nf not in range(2, 6 + 1):
+    if nf > 6:
         raise NotImplementedError("Selected nf is not implemented")
     nu = nf // 2
     return nu

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.
 
@@ -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
@@ -308,8 +310,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
@@ -322,7 +324,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])
@@ -521,7 +523,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
@@ -530,6 +532,8 @@ def compute_exact_alphaem_running(self, a_ref, nf, 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
@@ -576,12 +580,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 +660,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
@@ -668,6 +672,8 @@ def compute(self, a_ref, nf, 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
@@ -678,15 +684,15 @@ 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:
             # at the moment everything is expanded - and type has been checked in the constructor
             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 +704,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 +744,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 = matchings.lepton_number(seg.origin)
+                nlf = matchings.lepton_number(seg.target)
+                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(
+                        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

src/eko/evolution_operator/__init__.py

@@ -26,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__)
@@ -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, 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, 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, lepton_number(mu2_to), L, alphaem_running
             )
         ker = qed_ns.dispatcher(
             order,
@@ -704,40 +704,20 @@ 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
-                ]
+                [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.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, 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

src/eko/matchings.py

@@ -3,8 +3,10 @@
 from dataclasses import dataclass
 from typing import List, Union
 
+import numba as nb
 import numpy as np
 
+from .constants import MTAU
 from .io.types import EvolutionPoint as EPoint
 from .io.types import FlavorIndex, FlavorsNumber, SquaredScale
 from .quantities.heavy_quarks import MatchingScales
@@ -209,3 +211,23 @@ 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
+
+
+@nb.njit(cache=True)
+def lepton_number(q2):
+    """Compute the number of leptons.
+
+    Note: muons and electrons are always massless as for up, down and strange.
+
+    Parameters
+    ----------
+    q2 : float
+        scale
+
+    Returns
+    -------
+    int :
+       Number of leptons
+
+    """
+    return 3 if q2 > MTAU**2 else 2