X space conversion of the N3LO matching

extras/ome_n3lo/.gitignore

@@ -0,0 +1,3 @@
+*.pdf
+*.f
+*.txt

extras/ome_n3lo/convert_ome_xspace.py

@@ -0,0 +1,100 @@
+"""Dump a fast  x-space grid of the N3LO transition matrix elements.
+The output file have the structure: x_grid, nf=3, nf=4, nf=5.
+"""
+
+import numpy as np
+from click import progressbar
+from eko.mellin import Path
+from ekore.harmonics import cache as c
+from ekore.operator_matrix_elements.unpolarized.space_like import as3
+from eko.interpolation import lambertgrid
+from scipy import integrate
+
+from large_n_limit import Agg_asymptotic, Aqq_asymptotic
+
+XGRID = lambertgrid(500, 1e-6)
+"""X-grid."""
+
+LOG = 0
+"""Matching threshold displaced ?"""
+
+MAP_ENTRIES = {
+    "gg": (0, 0),
+    "gq": (0, 1),
+    "qg": (1, 0),
+    "qq": (1, 1),
+    "Hg": (2, 0),
+    "Hq": (2, 1),
+    "gH": (0, 2),
+    "HH": (2, 2),
+    "qq_ns": (0, 0),
+}
+
+
+def compute_ome(nf, n, is_singlet):
+    """Get the correct ome from eko."""
+    cache = c.reset()
+    if is_singlet:
+        return as3.A_singlet(n, cache, nf, L=0)
+    else:
+        return as3.A_ns(n, cache, nf, L=0)
+
+
+def compute_xspace_ome(entry, nf, x_grid=XGRID):
+    """Compute the x-space transition matrix element, returns A^3(x)."""
+    mellin_cut = 5e-2
+    is_singlet = "ns" not in entry
+
+    def integrand(u, x):
+        """Mellin inversion integrand."""
+        path = Path(u, np.log(x), is_singlet)
+        integrand = path.prefactor * x ** (-path.n) * path.jac
+        if integrand == 0.0:
+            return 0.0
+
+        # compute the N space ome
+        ome_n = compute_ome(nf, path.n, is_singlet)
+        idx1, idx2 = MAP_ENTRIES[entry]
+        ome_n = ome_n[idx1, idx2]
+        # subtract the large-N limit for diagonal terms (ie local and singular bits)
+        if entry in ["qq_ns", "qq"]:
+            ome_n -= Aqq_asymptotic(path.n, nf)
+        elif entry == "gg":
+            ome_n -= Agg_asymptotic(path.n, nf)
+
+        # recombine everything
+        return np.real(ome_n * integrand)
+
+    ome_x = []
+    print(f"Computing operator matrix element {entry} @ pto: 3, nf: {nf}")
+    # loop on xgrid
+    with progressbar(x_grid) as bar:
+        for x in bar:
+            if x == 1:
+                ome_x.append(0)
+                continue
+            res = integrate.quad(
+                lambda u: integrand(u, x),
+                0.5,
+                1.0 - mellin_cut,
+                epsabs=1e-12,
+                epsrel=1e-6,
+                limit=200,
+                full_output=1,
+            )[0]
+            ome_x.append(res)
+
+    return np.array(ome_x)
+
+
+def save_files(entry, ome_x, xgrid=XGRID):
+    """Write the space reuslt in a txt file."""
+    fname = f"x_space/A_{entry}.txt"
+    np.savetxt(fname, np.concatenate(([xgrid], np.array(ome_x))).T)
+
+
+if __name__ == "__main__":
+    # non diagonal temrms
+    for k in ["qq_ns", "gg", "gq", "qg", "qq", "Hg", "Hq"]:
+        result = [compute_xspace_ome(k, nf) for nf in [3, 4, 5]]
+        save_files(k, result)

extras/ome_n3lo/large_n_limit.py

@@ -0,0 +1,38 @@
+"""This file contains the large-N limit of the diagonal Matrix elements.
+
+The expansions are obtained using the notebook Agg_Aqq_largex_expansion.nb.
+
+We note that:
+    * the limit og :math:`A_{qq}` is the same for non-singlet like and singlet-like expansions.
+    I.e. the local and singular part are the same
+    * the :math:`A_{qq,ps}` temr is vanishing in the large-x limit, i.e. it's only regular.
+"""
+from ekore.harmonics import S1
+
+
+def Aqq_asymptotic(n, nf):
+    """The N3LO quark-to-quark transition matrix element large-N limit."""
+    return (
+        (20.362519064781296 - 3.4050138869326796 * nf) * S1(n)
+        - 51.033843609253296
+        + 3.1144841058729096 * nf
+    )
+
+
+def Agg_asymptotic(n, nf):
+    """The N3LO gluon-to-gluon transition matrix element large-N limit.
+    Follwing :cite:`Ablinger:2022wbb`:
+        * the fist part contains the limit of eq. 2.6 (except for :math:`a_{gg}^{(3)}`)
+        * the second part comes from eq. 4.6 and 4.7.
+    """
+    Agg_asy_incomplete = (
+        (-669.1554507291286 + 41.84286985333757 * nf) * S1(n)
+        - 565.4465327471261
+        + 28.65462637880661 * nf
+    )
+    agg_asy = (
+        - 49.5041510989361 * (-14.442649813264895 + nf) * S1(n)
+        + 619.2420126046355
+        - 17.52475977636971 * nf
+    )
+    return agg_asy + Agg_asy_incomplete