Prepare extension of DGLAP to QED

benchmarks/eko/benchmark_ad.py

@@ -3,8 +3,8 @@
 import numpy as np
 import pytest
 
+import eko.anomalous_dimensions.as2 as ad_as2
 import eko.anomalous_dimensions.harmonics as h
-import eko.anomalous_dimensions.nlo as ad_nlo
 from eko.constants import CA, CF, TR
 
 
@@ -136,8 +136,8 @@ def check_gamma_1_pegasus(N, NF):
     P1NSP = CF * ((CF - CA / 2.0) * PNPA + CA * PNSB + TR * NF * PNSC)
     P1NSM = CF * ((CF - CA / 2.0) * PNMA + CA * PNSB + TR * NF * PNSC)
 
-    np.testing.assert_allclose(ad_nlo.gamma_nsp_1(N, NF), -P1NSP)
-    np.testing.assert_allclose(ad_nlo.gamma_nsm_1(N, NF), -P1NSM)
+    np.testing.assert_allclose(ad_as2.gamma_nsp(N, NF), -P1NSP)
+    np.testing.assert_allclose(ad_as2.gamma_nsm(N, NF), -P1NSM)
 
     NS = N * N
     NT = NS * N
@@ -254,7 +254,7 @@ def check_gamma_1_pegasus(N, NF):
     P1Sgq = (CF * CF * PGQA + CF * CA * PGQB + TR * NF * CF * PGQC) * 4.0
     P1Sgg = (CA * CA * PGGA + TR * NF * (CA * PGGB + CF * PGGC)) * 4.0
 
-    gS1 = ad_nlo.gamma_singlet_1(N, NF)
+    gS1 = ad_as2.gamma_singlet(N, NF)
     np.testing.assert_allclose(gS1[0, 0], -P1Sqq)
     np.testing.assert_allclose(gS1[0, 1], -P1Sqg)
     np.testing.assert_allclose(gS1[1, 0], -P1Sgq)

benchmarks/eko/benchmark_evol_to_unity.py

@@ -2,6 +2,7 @@
 import numpy as np
 import pytest
 
+from eko import basis_rotation as br
 from eko.evolution_operator import Operator
 from eko.evolution_operator.grid import OperatorGrid
 from eko.interpolation import InterpolatorDispatcher
@@ -76,38 +77,39 @@ def test_operator_grid(
         o.compute()
         o_back.compute()
 
-        dim = o_back.op_members["NS_v"].value.shape
-        for k in ["NS_v", "NS_m", "NS_p"]:
+        dim = o_back.op_members[(br.non_singlet_pids_map["nsV"], 0)].value.shape
+        for k in ["nsV", "ns-", "ns+"]:
             np.testing.assert_allclose(
-                o.op_members[k].value @ o_back.op_members[k].value,
+                o.op_members[(br.non_singlet_pids_map[k], 0)].value
+                @ o_back.op_members[(br.non_singlet_pids_map[k], 0)].value,
                 np.eye(dim[0]),
                 atol=7e-2,
             )
         # qq
         np.testing.assert_allclose(
-            o_back.op_members["S_qq"].value @ o.op_members["S_qq"].value
-            + o_back.op_members["S_qg"].value @ o.op_members["S_gq"].value,
+            o_back.op_members[(100, 100)].value @ o.op_members[(100, 100)].value
+            + o_back.op_members[(100, 21)].value @ o.op_members[(21, 100)].value,
             np.eye(dim[0]),
             atol=7e-2,
         )
         # qg
         np.testing.assert_allclose(
-            o_back.op_members["S_qq"].value @ o.op_members["S_qg"].value
-            + o_back.op_members["S_qg"].value @ o.op_members["S_gg"].value,
+            o_back.op_members[(100, 100)].value @ o.op_members[(100, 21)].value
+            + o_back.op_members[(100, 21)].value @ o.op_members[(21, 21)].value,
             np.zeros(dim),
             atol=7e-4,
         )
         # gg
         np.testing.assert_allclose(
-            o_back.op_members["S_gg"].value @ o.op_members["S_gg"].value
-            + o_back.op_members["S_gq"].value @ o.op_members["S_qg"].value,
+            o_back.op_members[(21, 21)].value @ o.op_members[(21, 21)].value
+            + o_back.op_members[(21, 100)].value @ o.op_members[(100, 21)].value,
             np.eye(dim[0]),
             atol=9e-2,
         )
         # gq
         np.testing.assert_allclose(
-            o_back.op_members["S_gg"].value @ o.op_members["S_gq"].value
-            + o_back.op_members["S_gq"].value @ o.op_members["S_qq"].value,
+            o_back.op_members[(21, 21)].value @ o.op_members[(21, 100)].value
+            + o_back.op_members[(21, 100)].value @ o.op_members[(100, 100)].value,
             np.zeros(dim),
             atol=2e-3,
         )
@@ -124,7 +126,7 @@ def test_operator_grid(
     #     ome.compute( q2, L)
     #     ome_back.compute(q2, L)
 
-    #     dim = ome.ome_members["S_qq"].value.shape
+    #     dim = ome.ome_members[(100, 100)].value.shape
     #     ome_tensor = np.zeros((3,3,dim[0],dim[0]))
     #     ome_tensor_back = ome_tensor
     #     idx_dict = dict(zip(["g", "q", "H"],[0,1,2]))

doc/source/refs.bib

@@ -420,6 +420,17 @@ @article{Bertone:2013vaa
     year          = {2014}
 }
 
+@phdthesis{Carrazza:2015dea,
+    author = "Carrazza, Stefano",
+    title = "{Parton distribution functions with QED corrections}",
+    eprint = "1509.00209",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-ph",
+    doi = "10.13130/carrazza-stefano_phd2015-07-06",
+    school = "Milan U.",
+    year = "2015"
+}
+
 @article{Vermaseren:1997fq,
     author = "Vermaseren, J. A. M. and Larin, S. A. and van Ritbergen, T.",
     title = "{The four loop quark mass anomalous dimension and the invariant quark mass}",

src/eko/anomalous_dimensions/__init__.py

@@ -20,7 +20,8 @@
 import numba as nb
 import numpy as np
 
-from . import harmonics, lo, nlo, nnlo
+from .. import basis_rotation as br
+from . import as1, as2, as3, harmonics
 
 
 @nb.njit("Tuple((c16[:,:],c16,c16,c16[:,:],c16[:,:]))(c16[:,:])", cache=True)
@@ -52,9 +53,9 @@ def exp_singlet(gamma_S):
 
     See Also
     --------
-        eko.anomalous_dimensions.lo.gamma_singlet_0 : :math:`\gamma_{S}^{(0)}(N)`
-        eko.anomalous_dimensions.nlo.gamma_singlet_1 : :math:`\gamma_{S}^{(1)}(N)`
-        eko.anomalous_dimensions.nnlo.gamma_singlet_2 : :math:`\gamma_{S}^{(2)}(N)`
+        eko.anomalous_dimensions.as1.gamma_singlet : :math:`\gamma_{S}^{(0)}(N)`
+        eko.anomalous_dimensions.as2.gamma_singlet : :math:`\gamma_{S}^{(1)}(N)`
+        eko.anomalous_dimensions.as3.gamma_singlet : :math:`\gamma_{S}^{(2)}(N)`
     """
     # compute eigenvalues
     det = np.sqrt(
@@ -71,7 +72,7 @@ def exp_singlet(gamma_S):
     return exp, lambda_p, lambda_m, e_p, e_m
 
 
-@nb.njit("c16[:](u1,string,c16,u1)", cache=True)
+@nb.njit("c16[:](u1,u2,c16,u1)", cache=True)
 def gamma_ns(order, mode, n, nf):
     r"""
     Computes the tower of the non-singlet anomalous dimensions
@@ -80,7 +81,7 @@ def gamma_ns(order, mode, n, nf):
     ----------
         order : int
             perturbative order
-        mode : "m" | "p" | "v"
+        mode : 10201 | 10101 | 10200
             sector identifier
         n : complex
             Mellin variable
@@ -94,37 +95,39 @@ def gamma_ns(order, mode, n, nf):
 
     See Also
     --------
-        eko.anomalous_dimensions.lo.gamma_ns_0 : :math:`\gamma_{ns}^{(0)}(N)`
-        eko.anomalous_dimensions.nlo.gamma_nsp_1 : :math:`\gamma_{ns,+}^{(1)}(N)`
-        eko.anomalous_dimensions.nlo.gamma_nsm_1 : :math:`\gamma_{ns,-}^{(1)}(N)`
-        eko.anomalous_dimensions.nnlo.gamma_nsp_2 : :math:`\gamma_{ns,+}^{(2)}(N)`
-        eko.anomalous_dimensions.nnlo.gamma_nsm_2 : :math:`\gamma_{ns,-}^{(2)}(N)`
-        eko.anomalous_dimensions.nnlo.gamma_nsv_2 : :math:`\gamma_{ns,v}^{(2)}(N)`
+        eko.anomalous_dimensions.as1.gamma_ns : :math:`\gamma_{ns}^{(0)}(N)`
+        eko.anomalous_dimensions.as2.gamma_nsp : :math:`\gamma_{ns,+}^{(1)}(N)`
+        eko.anomalous_dimensions.as2.gamma_nsm : :math:`\gamma_{ns,-}^{(1)}(N)`
+        eko.anomalous_dimensions.as3.gamma_nsp : :math:`\gamma_{ns,+}^{(2)}(N)`
+        eko.anomalous_dimensions.as3.gamma_nsm : :math:`\gamma_{ns,-}^{(2)}(N)`
+        eko.anomalous_dimensions.as3.gamma_nsv : :math:`\gamma_{ns,v}^{(2)}(N)`
     """
     # cache the s-es
     sx = np.full(1, harmonics.harmonic_S1(n))
     # now combine
     gamma_ns = np.zeros(order + 1, np.complex_)
-    gamma_ns[0] = lo.gamma_ns_0(n, sx[0])
+    gamma_ns[0] = as1.gamma_ns(n, sx[0])
     # NLO and beyond
     if order >= 1:
         # TODO: pass the necessary harmonics to nlo gammas
-        if mode == "p":
-            gamma_ns_1 = nlo.gamma_nsp_1(n, nf)
+        if mode == 10101:
+            gamma_ns_1 = as2.gamma_nsp(n, nf)
         # To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here
-        elif mode in ["m", "v"]:
-            gamma_ns_1 = nlo.gamma_nsm_1(n, nf)
+        elif mode in [10201, 10200]:
+            gamma_ns_1 = as2.gamma_nsm(n, nf)
+        else:
+            raise NotImplementedError("Non-singlet sector is not implemented")
         gamma_ns[1] = gamma_ns_1
     # NNLO and beyond
     if order >= 2:
         sx = np.append(sx, harmonics.harmonic_S2(n))
         sx = np.append(sx, harmonics.harmonic_S3(n))
-        if mode == "p":
-            gamma_ns_2 = -nnlo.gamma_nsp_2(n, nf, sx)
-        elif mode == "m":
-            gamma_ns_2 = -nnlo.gamma_nsm_2(n, nf, sx)
-        elif mode == "v":
-            gamma_ns_2 = -nnlo.gamma_nsv_2(n, nf, sx)
+        if mode == 10101:
+            gamma_ns_2 = -as3.gamma_nsp(n, nf, sx)
+        elif mode == 10201:
+            gamma_ns_2 = -as3.gamma_nsm(n, nf, sx)
+        elif mode == 10200:
+            gamma_ns_2 = -as3.gamma_nsv(n, nf, sx)
         gamma_ns[2] = gamma_ns_2
     return gamma_ns
 
@@ -150,9 +153,9 @@ def gamma_singlet(order, n, nf):
 
     See Also
     --------
-        eko.anomalous_dimensions.lo.gamma_singlet_0 : :math:`\gamma_{S}^{(0)}(N)`
-        eko.anomalous_dimensions.nlo.gamma_singlet_1 : :math:`\gamma_{S}^{(1)}(N)`
-        eko.anomalous_dimensions.nnlo.gamma_singlet_2 : :math:`\gamma_{S}^{(2)}(N)`
+        eko.anomalous_dimensions.as1.gamma_singlet : :math:`\gamma_{S}^{(0)}(N)`
+        eko.anomalous_dimensions.as2.gamma_singlet : :math:`\gamma_{S}^{(1)}(N)`
+        eko.anomalous_dimensions.as3.gamma_singlet : :math:`\gamma_{S}^{(2)}(N)`
     """
     # cache the s-es
     sx = np.full(1, harmonics.harmonic_S1(n))
@@ -161,10 +164,10 @@ def gamma_singlet(order, n, nf):
         sx = np.append(sx, harmonics.harmonic_S3(n))
 
     gamma_singlet = np.zeros((order + 1, 2, 2), np.complex_)
-    gamma_singlet[0] = lo.gamma_singlet_0(n, sx[0], nf)
+    gamma_singlet[0] = as1.gamma_singlet(n, sx[0], nf)
     if order >= 1:
-        gamma_singlet[1] = nlo.gamma_singlet_1(n, nf)
+        gamma_singlet[1] = as2.gamma_singlet(n, nf)
     if order == 2:
         sx = np.append(sx, harmonics.harmonic_S4(n))
-        gamma_singlet[2] = -nnlo.gamma_singlet_2(n, nf, sx)
+        gamma_singlet[2] = -as3.gamma_singlet(n, nf, sx)
     return gamma_singlet

src/eko/anomalous_dimensions/aem1.py

@@ -0,0 +1,93 @@
+# -*- coding: utf-8 -*-
+import numba as nb
+
+from .. import constants
+from . import as1
+
+
+@nb.njit("c16(c16)", cache=True)
+def gamma_phq(N):
+    """
+    Computes the leading-order photon-quark anomalous dimension
+
+    Implements Eq. (2.5) of :cite:`Carrazza:2015dea`.
+
+    Parameters
+    ----------
+      N : complex
+        Mellin moment
+
+    Returns
+    -------
+      gamma_phq : complex
+        Leading-order photon-quark anomalous dimension :math:`\\gamma_{\\gamma q}^{(0)}(N)`
+    """
+
+    return as1.gamma_gq(N) / constants.CF
+
+
+@nb.njit("c16(c16,u1)", cache=True)
+def gamma_qph(N, nf):
+    """
+    Computes the leading-order quark-photon anomalous dimension
+
+    Implements Eq. (2.5) of :cite:`Carrazza:2015dea`.
+    But adding the :math:`N_C` and the :math:`2n_f` factors from :math:`\\theta` inside the
+    definition of :math:`\\gamma_{q \\gamma}^{(0)}(N)`.
+
+    Parameters
+    ----------
+      N : complex
+        Mellin moment
+      nf : int
+        Number of active flavors
+
+    Returns
+    -------
+      gamma_qph : complex
+        Leading-order quark-photon anomalous dimension :math:`\\gamma_{q \\gamma}^{(0)}(N)`
+    """
+    return as1.gamma_qg(N, nf) / constants.TR * constants.NC
+
+
+@nb.njit("c16(u1)", cache=True)
+def gamma_phph(nf):
+    """
+    Computes the leading-order photon-photon anomalous dimension
+
+    Implements Eq. (2.5) of :cite:`Carrazza:2015dea`.
+
+    Parameters
+    ----------
+      nf : int
+        Number of active flavors
+
+    Returns
+    -------
+      gamma_phph : complex
+        Leading-order phton-photon anomalous dimension :math:`\\gamma_{\\gamma \\gamma}^{(0)}(N)`
+    """
+
+    return 2 / 3 * constants.NC * 2 * nf
+
+
+@nb.njit("c16(c16,c16)", cache=True)
+def gamma_ns(N, s1):
+    """
+    Computes the leading-order non-singlet QED anomalous dimension.
+
+    Implements Eq. (2.5) of :cite:`Carrazza:2015dea`.
+
+    Parameters
+    ----------
+      N : complex
+        Mellin moment
+      s1 : complex
+        S1(N)
+
+    Returns
+    -------
+      gamma_ns : complex
+        Leading-order non-singlet QED anomalous dimension :math:`\\gamma_{ns}^{(0)}(N)`
+    """
+    return as1.gamma_ns(N, s1) / constants.CF