Merge pull request #245 from NNPDF/tlu_ome

Time like matching conditions

benchmarks/FF_LHAPDF_bench.py

@@ -104,7 +104,7 @@ def benchmark_lo(self, ff_index, Q0=10, mugrid=(100,)):
 
         self.run([theory_card], [operator_card], [FF_sets_lo[ff_index]])
 
-    def benchmark_nlo(self, ff_index, Q0=10, mugrid=(100,)):
+    def benchmark_nlo(self, ff_index, Q0=10, mugrid=(1.65,)):
         theory_card = {
             **base_theory,
             "PTO": 1,

doc/source/refs.bib

@@ -894,3 +894,16 @@ @article{Almasy:2011eq
     pages = "133--152",
     year = "2012"
 }
+
+@article{Cacciari:2005ry,
+    author = "Cacciari, Matteo and Nason, Paolo and Oleari, Carlo",
+    title = "{Crossing heavy-flavor thresholds in fragmentation functions}",
+    eprint = "hep-ph/0504192",
+    archivePrefix = "arXiv",
+    reportNumber = "BICOCCA-FT-05-8, LPTHE-05-09",
+    doi = "10.1088/1126-6708/2005/10/034",
+    journal = "JHEP",
+    volume = "10",
+    pages = "034",
+    year = "2005"
+}

doc/source/theory/TimeLike.rst

@@ -19,6 +19,8 @@ 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`.
+We also implement the time-like matching conditions up to |NLO| in :class:`~ekore.operator_matrix_elements.unpolarized.time_like` which
+are based on :cite:`Cacciari:2005ry`. The time-like matching conditions for |NNLO| are not known.
 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/operator_matrix_element.py

@@ -168,7 +168,7 @@ def quad_ker(
                 A = ome_ps.A_singlet(order, ker_base.n, sx, nf, L, is_msbar, sx_ns)
         else:
             if is_time_like:
-                A = ome_ut.A_singlet(order, ker_base.n, sx, nf, L, is_msbar, sx_ns)
+                A = ome_ut.A_singlet(order, ker_base.n, L)
             else:
                 A = ome_us.A_singlet(order, ker_base.n, sx, nf, L, is_msbar, sx_ns)
     else:
@@ -180,7 +180,7 @@ def quad_ker(
                 A = ome_ps.A_non_singlet(order, ker_base.n, sx, nf, L)
         else:
             if is_time_like:
-                A = ome_ut.A_non_singlet(order, ker_base.n, sx, nf, L)
+                A = ome_ut.A_non_singlet(order, ker_base.n, L)
             else:
                 A = ome_us.A_non_singlet(order, ker_base.n, sx, nf, L)
 

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

@@ -1,4 +1,4 @@
-"""The unpolarized time-like Altarelli-Parisi splitting kernels."""
+"""The unpolarized, time-like Altarelli-Parisi splitting kernels."""
 
 import numba as nb
 import numpy as np

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

@@ -1,4 +1,4 @@
-"""The unpolarized |LO| Altarelli-Parisi splitting kernels."""
+"""The unpolarized, time-like |LO| Altarelli-Parisi splitting kernels."""
 
 import numba as nb
 import numpy as np

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

@@ -1,4 +1,4 @@
-"""The unpolarized time-like |NLO| Altarelli-Parisi splitting kernels."""
+"""The unpolarized, time-like |NLO| Altarelli-Parisi splitting kernels."""
 
 import numba as nb
 import numpy as np

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

@@ -1,4 +1,4 @@
-"""The unpolarized time-like |NNLO| Altarelli-Parisi splitting kernels."""
+"""The unpolarized, time-like |NNLO| Altarelli-Parisi splitting kernels."""
 
 import numba as nb
 import numpy as np

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

@@ -1,15 +1,60 @@
-r"""The polarized, time-like |OME|."""
+r"""The unpolarized, time-like |OME|."""
 
 import numba as nb
+import numpy as np
+
+from . import as1
 
 
 @nb.njit(cache=True)
-def A_non_singlet(_matching_order, _n, _sx, _nf, _L):
-    """Compute the non-singlet |OME|."""
-    raise NotImplementedError("Time-like is not yet implemented")
+def A_non_singlet(matching_order, N, L):
+    r"""Compute the non-singlet |OME|.
+
+    Parameters
+    ----------
+    matching_order : tuple(int, int)
+        perturbative matching order
+    N : complex
+        Mellin moment
+    L : float
+        :math:`\ln(\mu_F^2 / m_h^2)`
+
+    Returns
+    -------
+    A_non_singlet : numpy.ndarray
+        non-singlet |OME|
+
+    """
+    A_ns = np.zeros((matching_order[0], 2, 2), np.complex_)
+    if matching_order[0] >= 1:
+        A_ns[0] = as1.A_ns()
+    if matching_order[0] >= 2:
+        A_ns[1] = as1.A_ns()
+    return A_ns
 
 
 @nb.njit(cache=True)
-def A_singlet(_matching_order, _n, _sx, _nf, _L, _is_msbar, _sx_ns=None):
-    """Compute the singlet |OME|."""
-    raise NotImplementedError("Time-like is not yet implemented")
+def A_singlet(matching_order, N, L):
+    r"""Compute the singlet |OME|.
+
+    Parameters
+    ----------
+    matching_order : tuple(int, int)
+        perturbative matching order
+    N : complex
+        Mellin moment
+    L : float
+        :math:`\ln(\mu_F^2 / m_h^2)`
+
+    Returns
+    -------
+    A_singlet : numpy.ndarray
+        singlet |OME|
+
+    """
+    A_singlet = np.zeros((matching_order[0], 3, 3), np.complex_)
+    if matching_order[0] >= 1:
+        A_singlet[0] = as1.A_singlet(N, L)
+    if matching_order[0] >= 2:
+        A_singlet[1] = as1.A_singlet(N, L)
+    return A_singlet

src/ekore/operator_matrix_elements/unpolarized/time_like/as1.py

@@ -0,0 +1,101 @@
+"""The unpolarized, time-like |NLO| matching conditions."""
+
+import numba as nb
+import numpy as np
+
+from eko.constants import CF
+
+
+@nb.njit(cache=True)
+def A_hg(N, L):
+    r"""Compute the |NLO| heavy-gluon |OME|.
+
+    Implements :eqref:`27` from :cite:`Cacciari:2005ry`.
+
+    Parameters
+    ----------
+    N : complex
+        Mellin moment
+    L : float
+        :math:`\ln(\mu_F^2 / m_h^2)`
+
+    Returns
+    -------
+    A_hg : complex
+        |NLO| heavy-gluon |OME|
+        :math:`A_{hg}^{S,(1)}`
+
+    """
+    result = (
+        2
+        * CF
+        * (
+            (2 + N + N**2) / (N * (N**2 - 1)) * (L - 1)
+            + 4 / (N - 1) ** 2
+            - 4 / N**2
+            + 2 / (N + 1) ** 2
+        )
+    )
+    return result
+
+
+@nb.njit(cache=True)
+def A_gg(L):
+    r"""Compute the |NLO| gluon-gluon |OME|.
+
+    Implements the Mellin transform of
+    :eqref:`24` from :cite:`Cacciari:2005ry`.
+    It is identical to the one in
+    :class:`~ekore.operator_matrix_elements.unpolarized.space_like.as1`.
+
+    Parameters
+    ----------
+    L : float
+        :math:`\ln(\mu_F^2 / m_h^2)`
+
+    Returns
+    -------
+    A_gg : complex
+        |NLO| gluon-gluon |OME| :math:`A_{gg,H}^{S,(1)}`
+
+    """
+    return -2.0 / 3.0 * L
+
+
+@nb.njit(cache=True)
+def A_singlet(N, L):
+    r"""Compute the |NLO| singlet |OME|.
+
+    Parameters
+    ----------
+    N : complex
+        Mellin moment
+    L : float
+        :math:`\ln(\mu_F^2 / m_h^2)`
+
+    Returns
+    -------
+    A_singlet : numpy.ndarray
+        |NLO| singlet |OME|
+        :math:`A^{S,(1)}`
+
+    """
+    result = np.array(
+        [[A_gg(L), 0 + 0j, 0], [0 + 0j, 0, 0], [A_hg(N, L), 0, 0]], np.complex_
+    )
+    return result
+
+
+@nb.njit(cache=True)
+def A_ns():
+    r"""Compute the |NLO| non-singlet |OME|.
+
+    Returns
+    -------
+    A_ns : numpy.ndarray
+        |NLO| non-singlet |OME|
+        :math:`A^{S,(1)}`
+
+    """
+    result = np.array([[0, 0], [0, 0]], np.complex_)
+    return result

tests/eko/evolution_operator/test_ome.py

@@ -81,7 +81,7 @@ def test_build_ome_nlo():
 
 
 def test_quad_ker_errors():
-    for p, t in [(True, False), (False, True), (True, True)]:
+    for p, t in [(True, False), (True, True)]:
         for mode0, mode1 in [
             (21, br.matching_hplus_pid),
             (200, br.matching_hminus_pid),
@@ -128,63 +128,64 @@ def test_quad_ker(monkeypatch):
         lambda *args: np.array([zeros, zeros, zeros]),
     )
     for is_log in [True, False]:
-        res_ns = quad_ker(
-            u=0,
-            order=(3, 0),
-            mode0=200,
-            mode1=200,
-            is_log=is_log,
-            logx=0.123,
-            areas=np.zeros(3),
-            backward_method=None,
-            a_s=0.0,
-            nf=3,
-            L=0.0,
-            sv_mode=sv.Modes.expanded,
-            Lsv=0.0,
-            is_msbar=False,
-            is_polarized=False,
-            is_time_like=False,
-        )
-        np.testing.assert_allclose(res_ns, 1.0)
-        res_s = quad_ker(
-            u=0,
-            order=(3, 0),
-            mode0=100,
-            mode1=100,
-            is_log=is_log,
-            logx=0.123,
-            areas=np.zeros(3),
-            backward_method=None,
-            a_s=0.0,
-            nf=3,
-            L=0.0,
-            sv_mode=sv.Modes.expanded,
-            Lsv=0.0,
-            is_msbar=False,
-            is_polarized=False,
-            is_time_like=False,
-        )
-        np.testing.assert_allclose(res_s, 1.0)
-        res_s = quad_ker(
-            u=0,
-            order=(3, 0),
-            mode0=100,
-            mode1=21,
-            is_log=is_log,
-            logx=0.0,
-            areas=np.zeros(3),
-            backward_method=None,
-            a_s=0.0,
-            nf=3,
-            L=0.0,
-            sv_mode=sv.Modes.expanded,
-            Lsv=0.0,
-            is_msbar=False,
-            is_polarized=False,
-            is_time_like=False,
-        )
-        np.testing.assert_allclose(res_s, 0.0)
+        for order, p, t in [((3, 0), False, False), ((1, 0), False, True)]:
+            res_ns = quad_ker(
+                u=0,
+                order=order,
+                mode0=200,
+                mode1=200,
+                is_log=is_log,
+                logx=0.123,
+                areas=np.zeros(3),
+                backward_method=None,
+                a_s=0.0,
+                nf=3,
+                L=0.0,
+                sv_mode=sv.Modes.expanded,
+                Lsv=0.0,
+                is_msbar=False,
+                is_polarized=p,
+                is_time_like=t,
+            )
+            np.testing.assert_allclose(res_ns, 1.0)
+            res_s = quad_ker(
+                u=0,
+                order=order,
+                mode0=100,
+                mode1=100,
+                is_log=is_log,
+                logx=0.123,
+                areas=np.zeros(3),
+                backward_method=None,
+                a_s=0.0,
+                nf=3,
+                L=0.0,
+                sv_mode=sv.Modes.expanded,
+                Lsv=0.0,
+                is_msbar=False,
+                is_polarized=p,
+                is_time_like=t,
+            )
+            np.testing.assert_allclose(res_s, 1.0)
+            res_s = quad_ker(
+                u=0,
+                order=order,
+                mode0=100,
+                mode1=21,
+                is_log=is_log,
+                logx=0.0,
+                areas=np.zeros(3),
+                backward_method=None,
+                a_s=0.0,
+                nf=3,
+                L=0.0,
+                sv_mode=sv.Modes.expanded,
+                Lsv=0.0,
+                is_msbar=False,
+                is_polarized=p,
+                is_time_like=t,
+            )
+            np.testing.assert_allclose(res_s, 0.0)
 
     # test expanded intrisic inverse kernels
     labels = [(200, 200), *br.singlet_labels]