Merge pull request #414 from NNPDF/improve_pgg_n3lo

Add newer parametrization for P_gg N3LO and use `use_fhmruvv=True` as default

doc/source/refs.bib

@@ -1049,6 +1049,17 @@ @article{Moch:2023tdj
     year = "2024"
 }
 
+@article{Falcioni:2024qpd,
+    author = "Falcioni, G. and Herzog, F. and Moch, S. and Pelloni, A. and Vogt, A.",
+    title = "{Four-loop splitting functions in QCD -- The gluon-gluon case --}",
+    eprint = "2410.08089",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-ph",
+    reportNumber = "ZU-TH 47/24, DESY-24-144, LTH 1384",
+    month = "10",
+    year = "2024"
+}
+
 @article{Ablinger:2024xtt,
     author = {Ablinger, J. and Behring, A. and Bl\"umlein, J. and De Freitas, A. and von Manteuffel, A. and Schneider, C. and Sch\"onwald, K.},
     title = "{The non-first-order-factorizable contributions to the three-loop single-mass operator matrix elements $A_{Qg}^{(3)}$ and $\Delta A_{Qg}^{(3)}$}",

doc/source/theory/N3LO_ad.rst

@@ -7,6 +7,28 @@ Moreover, the analytical structure of these function is already known to be comp
 since in Mellin space they include harmonics sum up to weight 7, for which an
 analytical expression is not available.
 
+We provide two different types of approximations, depending on the key ``use_fhmruvv``.
+The theory card parameter ``n3lo_ad_variation=(gg, gq, qg, qq, nsp, nsm, nsv)``, 
+which is set to ``(0,0,0,0,0,0,0)`` as default, can be varied to obtain parametrization uncertainties.
+
+In particular:
+
+* ``use_fhmruvv = True`` (default option) adopts the parametrizations as provided 
+  in :cite:`Moch:2017uml,Falcioni:2023luc,Falcioni:2023vqq,Falcioni:2024xyt,Falcioni:2024qpd`. 
+  For each of the 7 splitting functions, the approximation error can be obtained by varying 
+  the entries of ``n3lo_ad_variation`` in the range: ``0``, 
+  cental value (default), ``1`` down variation, ``2`` up variation.
+
+* ``use_fhmruvv = False`` adopts an in-house parametrization, constructed as described below.
+  In this case only variations of the singlet sector are available. Currently one 
+  can vary ``n3lo_ad_variation`` in the range: ``0-19`` for ``gg``, ``0-15`` for ``gg``,
+  ``0-15`` for ``qg``, ``0-6`` for ``qq``.
+  Note these approximations will be no longer updated and are now deprecated. 
+
+
+In house approximation
+----------------------
+
 Here, we describe the various assumptions and limits used in order to reconstruct a parametrization
 that can approximate their contribution.
 In particular we take advantage of some known physical constrains,
@@ -21,7 +43,7 @@ In any case |N3LO| |DGLAP| evolution at small-x, especially for singlet-like PDF
 until the splitting function resummation is available up to |NNLL|.
 
 Non-singlet sector
-------------------
+^^^^^^^^^^^^^^^^^^
 
 In the non-singlet sector we construct a parametrization for
 :math:`\gamma_{ns,-}^{(3)},\gamma_{ns,-}^{(3)},\gamma_{ns,s}^{(3)}` where:
@@ -147,7 +169,7 @@ In |EKO| they are implemented as follows:
             It is checked that this contribution is much more smaller than the values of :math:`B_4`.
 
 Singlet sector
---------------
+^^^^^^^^^^^^^^
 
 In the singlet sector we construct a parametrization for
 :math:`\gamma_{gg}^{(3)},\gamma_{gq}^{(3)},\gamma_{qg}^{(3)},\gamma_{qq}^{(3)}` where:

src/eko/io/runcards.py

@@ -49,7 +49,7 @@ class TheoryCard(DictLike):
     n3lo_ad_variation: N3LOAdVariation
     """|N3LO| anomalous dimension variation: ``(gg, gq, qg, qq, nsp, nsm,
     nsv)``."""
-    use_fhmruvv: Optional[bool] = False
+    use_fhmruvv: Optional[bool] = True
     """If True use the |FHMRUVV| |N3LO| anomalous dimensions."""
     matching_order: Optional[Order] = None
     """Matching conditions perturbative order tuple, ``(QCD, QED)``."""
@@ -213,7 +213,7 @@ def new_theory(self):
         new["n3lo_ad_variation"] = old.get("n3lo_ad_variation", (0, 0, 0, 0, 0, 0, 0))
         # here PTO: 0 means truly LO, no QED matching is available so far.
         new["matching_order"] = old.get("PTO_matching", [old["PTO"], 0])
-        new["use_fhmruvv"] = old.get("use_fhmruvv", False)
+        new["use_fhmruvv"] = old.get("use_fhmruvv", True)
         return TheoryCard.from_dict(new)
 
     @property

src/ekobox/cards.py

@@ -24,7 +24,7 @@
     xif=1.0,
     n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0),
     matching_order=[0, 0],
-    use_fhmruvv=False,
+    use_fhmruvv=True,
 )
 
 _operator = dict(

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

@@ -25,7 +25,7 @@
 
 
 @nb.njit(cache=True)
-def gamma_ns(order, mode, n, nf, n3lo_ad_variation, use_fhmruvv=False):
+def gamma_ns(order, mode, n, nf, n3lo_ad_variation, use_fhmruvv=True):
     r"""Compute the tower of the non-singlet anomalous dimensions.
 
     Parameters
@@ -100,7 +100,7 @@ def gamma_ns(order, mode, n, nf, n3lo_ad_variation, use_fhmruvv=False):
 
 
 @nb.njit(cache=True)
-def gamma_singlet(order, n, nf, n3lo_ad_variation, use_fhmruvv=False):
+def gamma_singlet(order, n, nf, n3lo_ad_variation, use_fhmruvv=True):
     r"""Compute the tower of the singlet anomalous dimensions matrices.
 
     Parameters
@@ -137,7 +137,7 @@ def gamma_singlet(order, n, nf, n3lo_ad_variation, use_fhmruvv=False):
 
 
 @nb.njit(cache=True)
-def gamma_ns_qed(order, mode, n, nf, n3lo_ad_variation, use_fhmruvv=False):
+def gamma_ns_qed(order, mode, n, nf, n3lo_ad_variation, use_fhmruvv=True):
     r"""Compute the grid of the QED non-singlet anomalous dimensions.
 
     Parameters
@@ -288,7 +288,7 @@ def choose_ns_ad_aem2(mode, n, nf, cache):
 
 
 @nb.njit(cache=True)
-def gamma_singlet_qed(order, n, nf, n3lo_ad_variation, use_fhmruvv=False):
+def gamma_singlet_qed(order, n, nf, n3lo_ad_variation, use_fhmruvv=True):
     r"""Compute the grid of the QED singlet anomalous dimensions matrices.
 
     Parameters
@@ -331,7 +331,7 @@ def gamma_singlet_qed(order, n, nf, n3lo_ad_variation, use_fhmruvv=False):
 
 
 @nb.njit(cache=True)
-def gamma_valence_qed(order, n, nf, n3lo_ad_variation, use_fhmruvv=False):
+def gamma_valence_qed(order, n, nf, n3lo_ad_variation, use_fhmruvv=True):
     r"""Compute the grid of the QED valence anomalous dimensions matrices.
 
     Parameters

src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/ggg.py

@@ -2,16 +2,30 @@
 :math:`\gamma_{gg}^{(3)}`."""
 
 import numba as nb
+import numpy as np
+
+from eko.constants import zeta3
 
 from ......harmonics import cache as c
-from ......harmonics.log_functions import lm11, lm12m1, lm13m1
+from ......harmonics.log_functions import (
+    lm11,
+    lm11m1,
+    lm11m2,
+    lm12m1,
+    lm12m2,
+    lm13m1,
+    lm13m2,
+    lm14m1,
+)
 
 
 @nb.njit(cache=True)
 def gamma_gg(n, nf, cache, variation):
     r"""Compute the |N3LO| gluon-gluon singlet anomalous dimension.
 
-    The routine is taken from :cite:`Moch:2023tdj`.
+    The routine is taken from :cite:`Falcioni:2024qpd`.
+    A previous version based only on the lowest 10 moments
+    was given in :cite:`Moch:2023tdj`.
 
     Parameters
     ----------
@@ -33,6 +47,9 @@ def gamma_gg(n, nf, cache, variation):
     S1 = c.get(c.S1, cache, n)
     S2 = c.get(c.S2, cache, n)
     S3 = c.get(c.S3, cache, n)
+    S4 = c.get(c.S4, cache, n)
+    Lm11m1 = lm11m1(n, S1)
+    Lm12m1 = lm12m1(n, S1, S2)
 
     nf2 = nf * nf
     nf3 = nf * nf2
@@ -41,80 +58,136 @@ def gamma_gg(n, nf, cache, variation):
     A4gluon = 40880.330 - 11714.246 * nf + 440.04876 * nf2 + 7.3627750 * nf3
     B4gluon = 68587.64 - 18143.983 * nf + 423.81135 * nf2 + 9.0672154 * 0.1 * nf3
 
+    # The coefficient of delta(1-x), also called the virtual anomalous
+    # dimension. nf^0 and nf^1 are still approximate, but the error at
+    # nf^1 is far too small to be relevant in this context.
+    if variation == 1:
+        B4gluon = B4gluon - 0.2
+    elif variation == 2:
+        B4gluon = B4gluon + 0.2
+
     Ccoeff = 8.5814120 * 10**4 - 1.3880515 * 10**4 * nf + 1.3511111 * 10**2 * nf2
     Dcoeff = 5.4482808 * 10**4 - 4.3411337 * 10**3 * nf - 2.1333333 * 10 * nf2
 
+    x1L4cff = 5.6460905 * 10 * nf - 3.6213992 * nf2
+    x1L3cff = 2.4755054 * 10**2 * nf - 4.0559671 * 10 * nf2 + 1.5802469 * nf3
+
     # The known coefficients of 1/x*ln^a x terms, a = 3,2
-    bfkl0 = -8.308617314 * 10**3
-    bfkl1 = -1.069119905 * 10**5 - 9.963830436 * 10**2 * nf
+    bfkl0 = -8.3086173 * 10**3
+    bfkl1 = -1.0691199 * 10**5 - 9.9638304 * 10**2 * nf
+
+    x0L6cff = 1.44 * 10**2 - 2.7786008 * 10 * nf + 7.9012346 * 0.1 * nf2
+    x0L5cff = -1.44 * 10**2 - 1.6208066 * 10**2 * nf + 1.4380247 * 10 * nf2
+    x0L4cff = (
+        2.6165784 * 10**4
+        - 3.3447551 * 10**3 * nf
+        + 9.1522635 * 10 * nf2
+        - 1.9753086 * 0.1 * nf3
+    )
 
     #  The resulting part of the function
     P3gg01 = (
         +bfkl0 * (-(6 / (-1 + n) ** 4))
         + bfkl1 * 2 / (-1 + n) ** 3
+        + x0L6cff * 720 / n**7
+        + x0L5cff * -120 / n**6
+        + x0L4cff * 24 / n**5
         + A4gluon * (-S1)
         + B4gluon
         + Ccoeff * lm11(n, S1)
         + Dcoeff * 1 / n
+        + x1L4cff * lm14m1(n, S1, S2, S3, S4)
+        + x1L3cff * lm13m1(n, S1, S2, S3)
     )
 
     # The selected approximations for nf = 3, 4, 5
     if nf == 3:
         P3ggApp1 = (
             P3gg01
-            + 3.4 * bfkl1 * -(1 / (-1 + n) ** 2)
-            - 345063.0 * 1 / ((-1 + n) * n)
-            + 86650.0 * (1 / n - 1 / (1 + n) + 1 / (2 + n) - 1 / (3 + n))
-            + 158160.0 * (-(1 / n**2))
-            - 15741.0 * lm12m1(n, S1, S2)
-            - 9417.0 * lm13m1(n, S1, S2, S3)
+            - 421311.0 * (-(1 / (-1 + n) ** 2) + 1 / n**2)
+            - 325557.0 * 1 / ((-1 + n) * n)
+            + 1679790.0 * (1 / (n + n**2))
+            - 1456863.0 * (1 / (2 + 3 * n + n**2))
+            + 3246307.0 * (-(1 / n**2) + 1 / (1 + n) ** 2)
+            + 2026324.0 * 2 / n**3
+            + 549188.0 * (-(6 / n**4))
+            + 8337.0 * Lm11m1
+            + 26718.0 * Lm12m1
+            - 27049.0 * lm13m2(n, S1, S2, S3)
         )
         P3ggApp2 = (
             P3gg01
-            + 5.4 * bfkl1 * -(1 / (-1 + n) ** 2)
-            - 1265632.0 * 1 / ((-1 + n) * n)
-            - 656644.0 * (1 / n - 1 / (1 + n) + 1 / (2 + n) - 1 / (3 + n))
-            - 1352233.0 * (-(1 / n**2))
-            + 203298.0 * lm12m1(n, S1, S2)
-            + 39112.0 * lm13m1(n, S1, S2, S3)
+            - 700113.0 * (-(1 / (-1 + n) ** 2) + 1 / n**2)
+            - 2300581.0 * 1 / ((-1 + n) * n)
+            + 896407.0 * (1 / n - n / (2 + 3 * n + n**2))
+            - 162733.0 * (1 / (6 + 5 * n + n**2))
+            - 2661862.0 * (-(1 / n**2) + 1 / (1 + n) ** 2)
+            + 196759.0 * 2 / n**3
+            - 260607.0 * (-(6 / n**4))
+            + 84068.0 * Lm11m1
+            + 346318.0 * Lm12m1
+            + 315725.0
+            * (
+                -3 * S1**2
+                + n * S1 * (np.pi**2 - 6 * S2)
+                - 3 * (S2 + 2 * n * (S3 - zeta3))
+            )
+            / (3 * n**2)
         )
     elif nf == 4:
         P3ggApp1 = (
             P3gg01
-            + 3.4 * bfkl1 * -(1 / (-1 + n) ** 2)
-            - 342625.0 * 1 / ((-1 + n) * n)
-            + 100372.0 * (1 / n - 1 / (1 + n) + 1 / (2 + n) - 1 / (3 + n))
-            + 189167.0 * (-(1 / n**2))
-            - 29762.0 * lm12m1(n, S1, S2)
-            - 12102.0 * lm13m1(n, S1, S2, S3)
+            - 437084.0 * (-(1 / (-1 + n) ** 2) + 1 / n**2)
+            - 361570.0 * 1 / ((-1 + n) * n)
+            + 1696070.0 * (1 / (n + n**2))
+            - 1457385.0 * (1 / (2 + 3 * n + n**2))
+            + 3195104.0 * (-(1 / n**2) + 1 / (1 + n) ** 2)
+            + 2009021.0 * 2 / n**3
+            + 544380.0 * (-(6 / n**4))
+            + 9938.0 * Lm11m1
+            + 24376.0 * Lm12m1
+            - 22143.0 * lm13m2(n, S1, S2, S3)
         )
         P3ggApp2 = (
             P3gg01
-            + 5.4 * bfkl1 * -(1 / (-1 + n) ** 2)
-            - 1271540.0 * 1 / ((-1 + n) * n)
-            - 649661.0 * (1 / n - 1 / (1 + n) + 1 / (2 + n) - 1 / (3 + n))
-            - 1334919.0 * (-(1 / n**2))
-            + 191263.0 * lm12m1(n, S1, S2)
-            + 36867.0 * lm13m1(n, S1, S2, S3)
+            - 706649.0 * (-(1 / (-1 + n) ** 2) + 1 / n**2)
+            - 2274637.0 * 1 / ((-1 + n) * n)
+            + 836544.0 * (1 / n - n / (2 + 3 * n + n**2))
+            - 199929.0 * (1 / (6 + 5 * n + n**2))
+            - 2683760.0 * (-(1 / n**2) + 1 / (1 + n) ** 2)
+            + 168802.0 * 2 / n**3
+            - 250799.0 * (-(6 / n**4))
+            + 36967.0 * Lm11m1
+            + 24530.0 * Lm12m1
+            - 71470.0 * lm12m2(n, S1, S2)
         )
     elif nf == 5:
         P3ggApp1 = (
             P3gg01
-            + 3.4 * bfkl1 * -(1 / (-1 + n) ** 2)
-            - 337540.0 * 1 / ((-1 + n) * n)
-            + 119366.0 * (1 / n - 1 / (1 + n) + 1 / (2 + n) - 1 / (3 + n))
-            + 223769.0 * (-(1 / n**2))
-            - 45129.0 * lm12m1(n, S1, S2)
-            - 15046.0 * lm13m1(n, S1, S2, S3)
+            - 439426.0 * (-(1 / (-1 + n) ** 2) + 1 / n**2)
+            - 293679.0 * 1 / ((-1 + n) * n)
+            + 1916281.0 * (1 / (n + n**2))
+            - 1615883.0 * (1 / (2 + 3 * n + n**2))
+            + 3648786.0 * (-(1 / n**2) + 1 / (1 + n) ** 2)
+            + 2166231.0 * 2 / n**3
+            + 594588.0 * (-(6 / n**4))
+            + 50406.0 * Lm11m1
+            + 24692.0 * Lm12m1
+            + 174067.0 * lm11m2(n, S1)
         )
         P3ggApp2 = (
             P3gg01
-            + 5.4 * bfkl1 * -(1 / (-1 + n) ** 2)
-            - 1274800.0 * 1 / ((-1 + n) * n)
-            - 637406.0 * (1 / n - 1 / (1 + n) + 1 / (2 + n) - 1 / (3 + n))
-            - 1314010.0 * (-(1 / n**2))
-            + 177882.0 * lm12m1(n, S1, S2)
-            + 34362.0 * lm13m1(n, S1, S2, S3)
+            - 705978.0 * (-(1 / (-1 + n) ** 2) + 1 / n**2)
+            - 2192234.0 * 1 / ((-1 + n) * n)
+            + 1730508.0 * (1 / (2 + 3 * n + n**2))
+            + 353143.0
+            * ((12 + 9 * n + n**2) / (6 * n + 11 * n**2 + 6 * n**3 + n**4))
+            - 2602682.0 * (-(1 / n**2) + 1 / (1 + n) ** 2)
+            + 178960.0 * 2 / n**3
+            - 218133.0 * (-(6 / n**4))
+            + 2285.0 * Lm11m1
+            + 19295.0 * Lm12m1
+            - 13719.0 * lm12m2(n, S1, S2)
         )
     else:
         raise NotImplementedError("nf=6 is not available at N3LO")

src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/ggq.py

@@ -23,7 +23,8 @@ def gamma_gq(n, nf, cache, variation):
     r"""Compute the |N3LO| gluon-quark singlet anomalous dimension.
 
     The routine is taken from :cite:`Falcioni:2024xyt`.
-    Lower moments were published also in :cite:`Moch:2023tdj`.
+    A previous version based only on the lowest 10 moments
+    was given in :cite:`Moch:2023tdj`.
 
 
     Parameters