Time Like Splitting Functions

benchmarks/FF_LHAPDF_bench.py

@@ -0,0 +1,149 @@
+"""Benchmark FFs from LHAPDF"""
+
+from banana import register
+
+from eko import interpolation
+from ekomark.benchmark.runner import Runner
+
+register(__file__)
+
+base_operator = {"ev_op_iterations": 10, "backward_inversion": "exact"}
+
+base_theory = {
+    "Qref": 91.1876,
+    "mc": 1.51,
+    "mb": 4.92,
+    "mt": 172.5,
+    "kcThr": 1.0,
+    "kbThr": 1.0,
+    "ktThr": 1.0,
+    "alphas": 0.118000,
+    "alphaqed": 0.007496,
+    "FNS": "ZM-VFNS",
+    "ModEv": "EXA",
+#    "ModEv": "TRN",
+}
+
+FF_sets_lo = [
+    "NNFF10_PIm_lo",
+    "NNFF10_PIp_lo",
+    "NNFF10_PIsum_lo",
+    "NNFF10_KAm_lo",
+    "NNFF10_KAp_lo",
+    "NNFF10_KAsum_lo",
+    "NNFF10_PRm_lo",
+    "NNFF10_PRp_lo",
+    "NNFF10_PRsum_lo",
+]
+FF_sets_nlo = [
+    "NNFF10_PIm_nlo",
+    "NNFF10_PIp_nlo",
+    "NNFF10_PIsum_nlo",
+    "NNFF10_KAm_nlo",
+    "NNFF10_KAp_nlo",
+    "NNFF10_KAsum_nlo",
+    "NNFF10_PRm_nlo",
+    "NNFF10_PRp_nlo",
+    "NNFF10_PRsum_nlo",
+    "MAPFF10NLOPIm",
+    "MAPFF10NLOPIp",
+    "MAPFF10NLOPIsum",
+    "MAPFF10NLOKAm",
+    "MAPFF10NLOKAp",
+    "MAPFF10NLOKAsum",
+]
+FF_sets_nnlo = [
+    "NNFF10_PIm_nnlo",
+    "NNFF10_PIp_nnlo",
+    "NNFF10_PIsum_nnlo",
+    "NNFF10_KAm_nnlo",
+    "NNFF10_KAp_nnlo",
+    "NNFF10_KAsum_nnlo",
+    "NNFF10_PRm_nnlo",
+    "NNFF10_PRp_nnlo",
+    "NNFF10_PRsum_nnlo",
+    "MAPFF10NNLOPIm",
+    "MAPFF10NNLOPIp",
+    "MAPFF10NNLOPIsum",
+    "MAPFF10NNLOKAm",
+    "MAPFF10NNLOKAp",
+    "MAPFF10NNLOKAsum",
+]
+
+
+class BenchmarkFF(Runner):
+    external = "LHAPDF"
+    rotate_to_evolution_basis = True
+
+    def skip_pdfs(self, _theory):
+        return [
+            -6,
+            6,
+            22,
+            "ph",
+            "T35",
+            "V35",
+        ]
+
+    def benchmark_lo(self, ff_index, Q0=10, mugrid=(100,)):
+        theory_card = {
+            **base_theory,
+            "PTO": 0,
+            "QED": 0,
+            "Q0": Q0,
+            "MaxNfPdf": 5,
+            "MaxNfAs": 5,
+        }
+
+        operator_card = {
+            **base_operator,
+            "mugrid": list(mugrid),
+            "time_like": True,
+            "interpolation_xgrid": interpolation.lambertgrid(100, 0.01),
+        }
+
+        self.run([theory_card], [operator_card], [FF_sets_lo[ff_index]])
+
+    def benchmark_nlo(self, ff_index, Q0=10, mugrid=(100,)):
+        theory_card = {
+            **base_theory,
+            "PTO": 1,
+            "QED": 0,
+            "Q0": Q0,
+            "MaxNfPdf": 5,
+            "MaxNfAs": 5,
+        }
+
+        operator_card = {
+            **base_operator,
+            "mugrid": list(mugrid),
+            "time_like": True,
+            "interpolation_xgrid": interpolation.lambertgrid(100, 0.01),
+        }
+
+        self.run([theory_card], [operator_card], [FF_sets_nlo[ff_index]])
+
+    def benchmark_nnlo(self, ff_index, Q0=10, mugrid=(100,)):
+        theory_card = {
+            **base_theory,
+            "PTO": 2,
+            "QED": 0,
+            "Q0": Q0,
+            "MaxNfPdf": 5,
+            "MaxNfAs": 5,
+        }
+
+        operator_card = {
+            **base_operator,
+            "mugrid": list(mugrid),
+            "time_like": True,
+            "interpolation_xgrid": interpolation.lambertgrid(100, 0.01),
+        }
+
+        self.run([theory_card], [operator_card], [FF_sets_nnlo[ff_index]])
+
+
+if __name__ == "__main__":
+#    BenchmarkFF().benchmark_lo(7)
+    BenchmarkFF().benchmark_nlo(10)
+#    BenchmarkFF().benchmark_nnlo(10)

doc/source/index.rst

@@ -39,6 +39,7 @@ EKO is ...
     theory/DGLAP
     theory/N3LO_ad
     theory/Matching
+    theory/TimeLike
     theory/MHOU
 
     zzz-refs

doc/source/refs.bib

@@ -827,4 +827,70 @@ @article{Moch:2014sna
     volume = "889",
     pages = "351--400",
     year = "2014"
-}
+}
+
+@article{Mitov:2006wy,
+    author = "Mitov, Alexander and Moch, Sven-Olaf",
+    title = "{QCD Corrections to Semi-Inclusive Hadron Production in Electron-Positron Annihilation at Two Loops}",
+    eprint = "hep-ph/0604160",
+    archivePrefix = "arXiv",
+    reportNumber = "DESY-06-043",
+    doi = "10.1016/j.nuclphysb.2006.05.018",
+    journal = "Nucl. Phys. B",
+    volume = "751",
+    pages = "18--52",
+    year = "2006"
+}
+
+@article{Gluck:1992zx,
+    author = "Gluck, M. and Reya, E. and Vogt, A.",
+    title = "{Parton fragmentation into photons beyond the leading order}",
+    reportNumber = "DO-TH-92-23",
+    doi = "10.1103/PhysRevD.51.1427",
+    journal = "Phys. Rev. D",
+    volume = "48",
+    pages = "116",
+    year = "1993",
+    note = "[Erratum: Phys.Rev.D 51, 1427 (1995)]"
+}
+
+@article{Mitov:2006ic,
+    author = "Mitov, A. and Moch, S. and Vogt, A.",
+    title = "{Next-to-Next-to-Leading Order Evolution of Non-Singlet Fragmentation Functions}",
+    eprint = "hep-ph/0604053",
+    archivePrefix = "arXiv",
+    reportNumber = "DESY-06-036, DCPT-06-40, IPPP-06-20",
+    doi = "10.1016/j.physletb.2006.05.005",
+    journal = "Phys. Lett. B",
+    volume = "638",
+    pages = "61--67",
+    year = "2006"
+}
+
+@article{Moch:2007tx,
+    author = "Moch, S. and Vogt, A.",
+    title = "{On third-order timelike splitting functions and top-mediated Higgs decay into hadrons}",
+    eprint = "0709.3899",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-ph",
+    reportNumber = "DESY-07-155, SFB-CPP-07-53, LTH-758",
+    doi = "10.1016/j.physletb.2007.10.069",
+    journal = "Phys. Lett. B",
+    volume = "659",
+    pages = "290--296",
+    year = "2008"
+}
+
+@article{Almasy:2011eq,
+    author = "Almasy, A. A. and Moch, S. and Vogt, A.",
+    title = "{On the Next-to-Next-to-Leading Order Evolution of Flavour-Singlet Fragmentation Functions}",
+    eprint = "1107.2263",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-ph",
+    reportNumber = "LTH-916, DESY-11-108, SFB-CPP-11-33, LPN11-34",
+    doi = "10.1016/j.nuclphysb.2011.08.028",
+    journal = "Nucl. Phys. B",
+    volume = "854",
+    pages = "133--152",
+    year = "2012"
+}

doc/source/shared/abbreviations.rst

@@ -4,6 +4,9 @@
 .. |PDF| replace::
    :abbr:`PDF (Parton Distribution Function(s))`
 
+.. |FF| replace::
+   :abbr:`FF (Fragmentation Function(s))`
+
 .. FNS
 
 .. |FNS| replace::

doc/source/theory/TimeLike.rst

@@ -0,0 +1,23 @@
+Time-Like Evolution
+===================
+
+Due to confinement in |QCD| we can not observe partons, such as quarks and gluons,
+directly in particle collider experiments.
+Instead, stable hadrons are detected which originate from parton interactions.
+
+The fragmentation functions (|FF|) encode the information
+on the probability for a hadron carrying a specified momentum fraction to 'fragment'
+from a given parton. These functions are non-perturbative and usually require a global |QCD|
+analysis of experimental data involving different processes for their reliable
+determination. This makes the |FF| similar to |PDF| as both rely
+on similar factorization theorems and, thus, on similar |RGE|.
+In practice, the relevant Feynman diagrams can indeed be related by a crossing
+symmetry which in turn means certain Mandelstam variables become for |FF|
+time-like instead of space-like.
+The relevant setting in the operator card is thus called ``time_like = True``.
+
+We implement the time-like |DGLAP| anomalous dimensions up to |NNLO| in :class:`~ekore.anomalous_dimensions.unpolarized.time_like`.
+The implementation for the |LO| and |NLO| splitting functions is based on :cite:`Mitov:2006wy, Gluck:1992zx` and the implementation for
+the |NNLO| splitting functions is based on :cite:`Mitov:2006ic, Moch:2007tx, Almasy:2011eq`.
+Supplying new anomalous dimensions and new matching conditions is the only change required for the eko program (e.g. the
+solution strategies are unaffected).

src/eko/evolution_operator/grid.py

@@ -107,6 +107,7 @@ def __init__(
             logger.warning("Evolution: In LO we use the exact solution always!")
 
         logger.info(dict(polarized=configs.polarized))
+        logger.info(dict(time_like=configs.time_like))
 
         self.config = config
         self.q2_grid = mu2grid

src/ekore/anomalous_dimensions/unpolarized/time_like/__init__.py

@@ -1,21 +1,74 @@
-r"""The unpolarized, time-like Altarelli-Parisi splitting kernels.
+"""The unpolarized time-like Altarelli-Parisi splitting kernels."""
 
-Normalization is given by
+import numba as nb
+import numpy as np
 
-.. math::
-    \mathbf{P}(x) = \sum\limits_{j=0} a_s^{j+1} \mathbf P^{(j)}(x)
+from . import as1, as2, as3
 
-with :math:`a_s = \frac{\alpha_S(\mu^2)}{4\pi}`.
-"""
 
-import numba as nb
+@nb.njit(cache=True)
+def gamma_ns(order, mode, n, nf):
+    r"""Compute the tower of the non-singlet anomalous dimensions.
 
+    Parameters
+    ----------
+    order : tuple(int,int)
+        perturbative orders
+    mode : 10201 | 10101 | 10200
+        sector identifier
+    n : complex
+        Mellin variable
+    nf : int
+        Number of active flavors
 
-@nb.njit(cache=True)
-def gamma_ns(_order, _mode, _n, _nf):
-    raise NotImplementedError("Polarised is not yet implemented")
+    Returns
+    -------
+    numpy.ndarray
+        non-singlet anomalous dimensions
+
+    """
+    gamma_ns = np.zeros(order[0], np.complex_)
+    gamma_ns[0] = as1.gamma_ns(n)
+    if order[0] >= 2:
+        if mode == 10101:
+            gamma_ns_1 = as2.gamma_nsp(n, nf)
+        elif mode in [10201, 10200]:
+            gamma_ns_1 = as2.gamma_nsm(n, nf)
+        gamma_ns[1] = gamma_ns_1
+    if order[0] >= 3:
+        if mode == 10101:
+            gamma_ns_2 = as3.gamma_nsp(n, nf)
+        elif mode == 10201:
+            gamma_ns_2 = as3.gamma_nsm(n, nf)
+        elif mode == 10200:
+            gamma_ns_2 = as3.gamma_nsv(n, nf)
+        gamma_ns[2] = gamma_ns_2
+    return gamma_ns
 
 
 @nb.njit(cache=True)
-def gamma_singlet(_order, _n, _nf):
-    raise NotImplementedError("Polarised is not yet implemented")
+def gamma_singlet(order, n, nf):
+    r"""Compute the tower of the singlet anomalous dimensions' matrices.
+
+    Parameters
+    ----------
+    order : tuple(int,int)
+        perturbative orders
+    n : complex
+        Mellin variable
+    nf : int
+        Number of active flavors
+
+    Returns
+    -------
+    numpy.ndarray
+        singlet anomalous dimensions matrices
+
+    """
+    gamma_s = np.zeros((order[0], 2, 2), np.complex_)
+    gamma_s[0] = as1.gamma_singlet(n, nf)
+    if order[0] >= 2:
+        gamma_s[1] = as2.gamma_singlet(n, nf)
+    if order[0] >= 3:
+        gamma_s[2] = as3.gamma_singlet(n, nf)
+    return gamma_s