LO polarised anomalous dimension

benchmarks/eko/benchmark_evol_to_unity.py

@@ -55,8 +55,9 @@ def test_operator_grid(
         q20 = 30
         q21 = 50
         nf = 4
-        o = Operator(g.config, g.managers, nf, q20, q21)
-        o_back = Operator(g.config, g.managers, nf, q21, q20)
+        p = False
+        o = Operator(g.config, g.managers, nf, p, q20, q21)
+        o_back = Operator(g.config, g.managers, nf, p, q21, q20)
         o.compute()
         o_back.compute()
 

benchmarks/performance/harmonics.py

@@ -17,12 +17,12 @@ def setup(self):
 
     def time_as1_sing(self):
         for n in self.ns:
-            ad.gamma_singlet(0, n, NF)
+            ad.gamma_singlet(0, n, NF, p=False)
 
     def time_as2_sing(self):
         for n in self.ns:
-            ad.gamma_singlet(1, n, NF)
+            ad.gamma_singlet(1, n, NF, p=False)
 
     def time_as3_sing(self):
         for n in self.ns:
-            ad.gamma_singlet(2, n, NF)
+            ad.gamma_singlet(2, n, NF, p=False)

src/eko/anomalous_dimensions/__init__.py

@@ -21,7 +21,7 @@
 
 from .. import basis_rotation as br
 from .. import harmonics
-from . import aem1, aem2, as1, as1aem1, as2, as3, as4
+from . import aem1, aem2, as1, as1aem1, as2, as3, as4, asp1
 
 
 @nb.njit(cache=True)
@@ -73,7 +73,7 @@ def exp_singlet(gamma_S):
 
 
 @nb.njit(cache=True)
-def gamma_ns(order, mode, n, nf):
+def gamma_ns(order, mode, n, nf, p):
     r"""Computes the tower of the non-singlet anomalous dimensions
 
     Parameters
@@ -86,6 +86,8 @@ def gamma_ns(order, mode, n, nf):
         Mellin variable
     nf : int
         Number of active flavors
+    p : Boolean
+        Polarised (True) or un-Polarised (False)
 
     Returns
     -------
@@ -120,9 +122,12 @@ def gamma_ns(order, mode, n, nf):
         sx = harmonics.sx(n, max_weight=order[0] + 1)
     # now combine
     gamma_ns = np.zeros(order[0], np.complex_)
-    gamma_ns[0] = as1.gamma_ns(n, sx[0])
-    # NLO and beyond
-    if order[0] >= 2:
+    if p == False:
+        gamma_ns[0] = as1.gamma_ns(n, sx[0])
+    elif p == True:
+        gamma_ns[0] = asp1.gamma_pns(n, sx[0])
+    # NLO and beyondd
+    if order[0] >= 2 and p == False:
         if mode == 10101:
             gamma_ns_1 = as2.gamma_nsp(n, nf, sx)
         # To fill the full valence vector in NNLO we need to add gamma_ns^1 explicitly here
@@ -131,29 +136,41 @@ def gamma_ns(order, mode, n, nf):
         else:
             raise NotImplementedError("Non-singlet sector is not implemented")
         gamma_ns[1] = gamma_ns_1
+
+    elif order[0] >= 2 and p == True:
+        raise Exception("Polarised case is not known at this order")
+
     # NNLO and beyond
-    if order[0] >= 3:
+    if order[0] >= 3 and p == False:
         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
+
+    elif order[0] >= 3 and p == True:
+        raise Exception("Polarised case is not known at this order")
+
     # N3LO
-    if order[0] >= 4:
+    if order[0] >= 4 and p == False:
         if mode == 10101:
             gamma_ns_3 = as4.gamma_nsp(n, nf, full_sx_cache)
         elif mode == 10201:
             gamma_ns_3 = as4.gamma_nsm(n, nf, full_sx_cache)
         elif mode == 10200:
             gamma_ns_3 = as4.gamma_nsv(n, nf, full_sx_cache)
         gamma_ns[3] = gamma_ns_3
+
+    elif order[0] >= 4 and p == True:
+        raise Exception("Polarised case is not known at this order")
+
     return gamma_ns
 
 
 @nb.njit(cache=True)
-def gamma_singlet(order, n, nf):
+def gamma_singlet(order, n, nf, p):
     r"""Computes the tower of the singlet anomalous dimensions matrices
 
     Parameters
@@ -164,7 +181,8 @@ def gamma_singlet(order, n, nf):
         Mellin variable
     nf : int
         Number of active flavors
-
+    p : Boolean
+        Polarised (True) or un-Polarised (False)
     Returns
     -------
     numpy.ndarray
@@ -196,11 +214,25 @@ def gamma_singlet(order, n, nf):
         sx = harmonics.sx(n, max_weight=order[0])
 
     gamma_s = np.zeros((order[0], 2, 2), np.complex_)
-    gamma_s[0] = as1.gamma_singlet(n, sx[0], nf)
-    if order[0] >= 2:
+
+    if p == False:
+        gamma_s[0] = as1.gamma_singlet(n, sx[0], nf)
+    elif p == True:
+        gamma_s[0] = asp1.gamma_psinglet(n, sx[0], nf)
+
+    if order[0] >= 2 and p == False:
         gamma_s[1] = as2.gamma_singlet(n, nf, sx)
-    if order[0] >= 3:
+    elif order[0] >= 2 and p == True:
+        raise Exception("Polarised case is not known at this order")
+
+    if order[0] >= 3 and p == False:
         gamma_s[2] = as3.gamma_singlet(n, nf, sx)
-    if order[0] >= 4:
+    elif order[0] >= 3 and p == True:
+        raise Exception("Polarised case is not known at this order")
+
+    if order[0] >= 4 and p == False:
         gamma_s[3] = as4.gamma_singlet(n, nf, full_sx_cache)
+    elif order[0] >= 4 and p == True:
+        raise Exception("Polarised case is not known at this order")
+
     return gamma_s

src/eko/anomalous_dimensions/aem1.py

@@ -8,6 +8,7 @@
 from . import as1
 
 
+
 @nb.njit(cache=True)
 def gamma_phq(N):
     """

src/eko/anomalous_dimensions/asp1.py

@@ -0,0 +1,127 @@
+# -*- coding: utf-8 -*-
+"""This file contains the leading-order polarised  Altarelli-Parisi splitting kernels."""
+
+import numba as nb
+import numpy as np
+
+from .. import constants
+from .as1 import gamma_ns as gamma_pns
+
+# @nb.njit(cache=True)
+# def gamma_pns(N, s1):
+# Computes the leading-order non-singlet anomalous dimension for the polarised case.
+# This is going to be the same expression as the one for the unpolarised case.
+
+# gamma = -(3.0 - 4.0 * s1 + 2.0 / N / (N + 1.0))
+# result = constants.CF * gamma
+# return result
+
+
+@nb.njit(cache=True)
+def gamma_pqg(N, nf):
+    """
+    Computes the leading-order polarised quark-gluon anomalous dimension
+    #requires citation
+
+    Parameters
+    ----------
+      N : complex
+        Mellin moment
+      nf : int
+        Number of active flavors
+
+    Returns
+    -------
+      gamma_qg : complex
+        Leading-order polarised quark-gluon anomalous dimension :math:`\\gamma_{qg}^{(0)}(N)`
+    """
+    gamma = -(N - 1) / N / (N + 1)
+    result = 2.0 * constants.TR * 2.0 * nf * gamma
+    return result
+
+
+@nb.njit(cache=True)
+def gamma_pgq(N):
+    """
+    Computes the leading-order polarised gluon-quark anomalous dimension
+
+
+    Parameters
+    ----------
+      N : complex
+        Mellin moment
+
+    Returns
+    -------
+      gamma_gq : complex
+        Leading-order gluon-quark anomalous dimension :math:`\\gamma_{gq}^{(0)}(N)`
+    """
+    gamma = -(N + 2) / N / (N + 1)
+    result = 2.0 * constants.CF * gamma
+    return result
+
+
+@nb.njit(cache=True)
+def gamma_pgg(N, s1, nf):
+    """
+    Computes the leading-order polarised gluon-gluon anomalous dimension
+
+
+    Parameters
+    ----------
+      N : complex
+        Mellin moment
+      s1 : complex
+        harmonic sum :math:`S_{1}`
+      nf : int
+        Number of active flavors
+
+    Returns
+    -------
+      gamma_gg : complex
+        Leading-order gluon-gluon anomalous dimension :math:`\\gamma_{gg}^{(0)}(N)`
+    """
+    gamma = s1 - 2 / N / (N + 1)
+    result = constants.CA * (4.0 * gamma - 11.0 / 3.0) + 4.0 / 3.0 * constants.TR * nf
+    return result
+
+
+
+@nb.njit(cache=True)
+def gamma_psinglet(N, s1, nf):
+    r"""
+      Computes the leading-order polarised singlet anomalous dimension matrix
+
+      .. math::
+          \gamma_S^{(0)} = \left(\begin{array}{cc}
+            \gamma_{qq}^{(0)} & \gamma_{qg}^{(0)}\\
+            \gamma_{gq}^{(0)} & \gamma_{gg}^{(0)}
+          \end{array}\right)
+
+      Parameters
+      ----------
+        N : complex
+          Mellin moment
+        s1 : complex
+          harmonic sum :math:`S_{1}`
+        nf : int
+          Number of active flavors
+
+      Returns
+      -------
+        gamma_S_0 : numpy.ndarray
+          Leading-order singlet anomalous dimension matrix :math:`\gamma_{S}^{(0)}(N)`
+
+      See Also
+      --------
+        gamma_ns : :math:`\gamma_{qq}^{(0)}`
+        gamma_qg : :math:`\gamma_{qg}^{(0)}`
+        gamma_gq : :math:`\gamma_{gq}^{(0)}`
+        gamma_gg : :math:`\gamma_{gg}^{(0)}`
+    """
+    gamma_pqq = gamma_pns(N, s1)
+    gamma_pS_0 = np.array(
+        [[gamma_pqq, gamma_pqg(N, nf)], [gamma_pgq(N), gamma_pgg(N, s1, nf)]],
+        np.complex_,
+    )
+    return gamma_pS_0

src/eko/evolution_operator/__init__.py

@@ -11,15 +11,14 @@
 
 import numba as nb
 import numpy as np
-from scipy import integrate
-
 from .. import anomalous_dimensions as ad
 from .. import basis_rotation as br
 from .. import interpolation, mellin
 from .. import scale_variations as sv
 from ..kernels import non_singlet as ns
 from ..kernels import singlet as s
 from ..member import OpMember
+from scipy import integrate
 
 logger = logging.getLogger(__name__)
 
@@ -84,7 +83,7 @@ def path(self):
 
     @property
     def n(self):
-        """Returs the Mellin moment :math:`N`."""
+        """Return the Mellin moment :math:`N`."""
         return self.path.n
 
     def integrand(
@@ -125,6 +124,7 @@ def quad_ker(
     as0,
     as_raw,
     nf,
+    p,
     L,
     ev_op_iterations,
     ev_op_max_order,
@@ -159,6 +159,8 @@ def quad_ker(
         coupling value at the process scale
     nf : int
         number of active flavors
+    p : Boolean
+        Polarised (True) or un-Polarised (False)
     L : float
         logarithm of the squared ratio of factorization and renormalization scale
     ev_op_iterations : int
@@ -182,7 +184,7 @@ def quad_ker(
 
     # compute the actual evolution kernel
     if ker_base.is_singlet:
-        gamma_singlet = ad.gamma_singlet(order, ker_base.n, nf)
+        gamma_singlet = ad.gamma_singlet(order, ker_base.n, nf, p)
         # scale var exponentiated is directly applied on gamma
         if sv_mode == sv.Modes.exponentiated:
             gamma_singlet = sv.exponentiated.gamma_variation(
@@ -198,14 +200,14 @@ def quad_ker(
             ev_op_iterations,
             ev_op_max_order,
         )
-        # scale var expanded is applied on the kernel
+
         if sv_mode == sv.Modes.expanded and not is_threshold:
             ker = np.ascontiguousarray(
                 sv.expanded.singlet_variation(gamma_singlet, as_raw, order, nf, L)
             ) @ np.ascontiguousarray(ker)
         ker = select_singlet_element(ker, mode0, mode1)
     else:
-        gamma_ns = ad.gamma_ns(order, mode0, ker_base.n, nf)
+        gamma_ns = ad.gamma_ns(order, mode0, ker_base.n, nf, p)
         if sv_mode == sv.Modes.exponentiated:
             gamma_ns = sv.exponentiated.gamma_variation(gamma_ns, order, nf, L)
         ker = ns.dispatcher(
@@ -254,7 +256,14 @@ class Operator(sv.ModeMixin):
     full_labels = br.full_labels
 
     def __init__(
-        self, config, managers, nf, q2_from, q2_to, mellin_cut=5e-2, is_threshold=False
+        self,
+        config,
+        managers,
+        nf,
+        q2_from,
+        q2_to,
+        mellin_cut=5e-2,
+        is_threshold=False,
     ):
         self.config = config
         self.managers = managers
@@ -377,7 +386,8 @@ def quad_ker(self, label, logx, areas):
             as0=self.a_s[0],
             as_raw=self.a_s[2],
             nf=self.nf,
-            L=np.log(self.xif2),
+            p=self.config["p"],
+            L=np.log(self.fact_to_ren),
             ev_op_iterations=self.config["ev_op_iterations"],
             ev_op_max_order=tuple(self.config["ev_op_max_order"]),
             sv_mode=self.sv_mode,