diff --git a/doc/source/refs.bib b/doc/source/refs.bib index 769f2e326..7c820ec93 100644 --- a/doc/source/refs.bib +++ b/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)}$}", diff --git a/doc/source/theory/N3LO_ad.rst b/doc/source/theory/N3LO_ad.rst index 1fb13a15e..3a9141203 100644 --- a/doc/source/theory/N3LO_ad.rst +++ b/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: diff --git a/src/eko/io/runcards.py b/src/eko/io/runcards.py index 9de351c19..b83cf0db8 100644 --- a/src/eko/io/runcards.py +++ b/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 diff --git a/src/ekobox/cards.py b/src/ekobox/cards.py index 295993634..7e9774c7a 100644 --- a/src/ekobox/cards.py +++ b/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( diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py index cf569d9f7..167110115 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/__init__.py +++ b/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 diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/ggg.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/ggg.py index 230e94f11..71b29d988 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/ggg.py +++ b/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") diff --git a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/ggq.py b/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/ggq.py index 1d603a1ee..1affbff95 100644 --- a/src/ekore/anomalous_dimensions/unpolarized/space_like/as4/fhmruvv/ggq.py +++ b/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 diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py index ed003f925..69cc29e5a 100644 --- a/tests/eko/evolution_operator/test_init.py +++ b/tests/eko/evolution_operator/test_init.py @@ -45,7 +45,7 @@ def test_quad_ker_errors(): is_polarized=True, is_time_like=True, n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0), - use_fhmruvv=False, + use_fhmruvv=True, ) @@ -112,7 +112,7 @@ def test_quad_ker(monkeypatch): is_polarized=p, is_time_like=t, n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0), - use_fhmruvv=False, + use_fhmruvv=True, ) np.testing.assert_allclose(res_ns, res) for label in [(br.non_singlet_pids_map["ns+"], 0), (100, 100)]: @@ -141,7 +141,7 @@ def test_quad_ker(monkeypatch): is_polarized=polarized, is_time_like=False, n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0), - use_fhmruvv=False, + use_fhmruvv=True, ) np.testing.assert_allclose(res_sv, 1.0) for label in [ @@ -177,7 +177,7 @@ def test_quad_ker(monkeypatch): n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0), is_polarized=False, is_time_like=False, - use_fhmruvv=False, + use_fhmruvv=True, ) np.testing.assert_allclose(res_sv, 1.0) @@ -205,7 +205,7 @@ def test_quad_ker(monkeypatch): n3lo_ad_variation=(0, 0, 0, 0, 0, 0, 0), is_polarized=False, is_time_like=False, - use_fhmruvv=False, + use_fhmruvv=True, ) np.testing.assert_allclose(res_ns, 0.0) diff --git a/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as4_fhmv.py b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as4_fhmv.py index a6671e3b2..4f778c0dd 100644 --- a/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as4_fhmv.py +++ b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_as4_fhmv.py @@ -50,8 +50,8 @@ def test_momentum_conservation(): ) np.testing.assert_allclose( g_singlet[:, 0, 1] + g_singlet[:, 1, 1], - [-0.134766, 0.465174, -0.734706], - atol=2e-5, + [-0.0300842, +0.283004, -0.343172], + atol=1e-6, ) diff --git a/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py index 5fdc1ee06..f726f5248 100644 --- a/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py +++ b/tests/ekore/anomalous_dimensions/unpolarized/space_like/test_init.py @@ -118,7 +118,7 @@ def test_gamma_ns(): # as4 assert_allclose( ad_us.gamma_ns( - (4, 0), br.non_singlet_pids_map["ns-"], 1, nf, n3lo_ad_variation + (4, 0), br.non_singlet_pids_map["ns-"], 1, nf, n3lo_ad_variation, use_fhmruvv = False ), np.zeros(4), atol=2e-4, @@ -126,7 +126,7 @@ def test_gamma_ns(): # N3LO valence has a spurious pole, need to add a small shift assert_allclose( ad_us.gamma_ns( - (4, 0), br.non_singlet_pids_map["nsV"], 1 + 1e-6, nf, n3lo_ad_variation + (4, 0), br.non_singlet_pids_map["nsV"], 1 + 1e-6, nf, n3lo_ad_variation, use_fhmruvv = False ), np.zeros(4), atol=5e-4,