N3LO ad variations

benchmarks/lha_paper_bench.py

@@ -248,8 +248,7 @@ def run_lha(self, theory_updates):
     # Benchmark to LHA
     obj = BenchmarkFFNS_polarized()
     # obj = BenchmarkFFNS()
-
-    # obj.benchmark_plain(0)
+    # obj.benchmark_plain(1)
     obj.benchmark_sv(1)
 
     # # VFNS benchmarks with LHA settings

doc/source/refs.bib

@@ -830,3 +830,54 @@ @article{Falcioni:2023luc
     month = "2",
     year = "2023"
 }
+@article{Gluck:1995yr,
+    author = "Gluck, M. and Reya, E. and Stratmann, M. and Vogelsang, W.",
+    title = "{Next-to-leading order radiative parton model analysis of polarized deep inelastic lepton - nucleon scattering}",
+    eprint = "hep-ph/9508347",
+    archivePrefix = "arXiv",
+    reportNumber = "DO-TH-95-13, RAL-TR-95-042",
+    doi = "10.1103/PhysRevD.53.4775",
+    journal = "Phys. Rev. D",
+    volume = "53",
+    pages = "4775--4786",
+    year = "1996"
+}
+
+@article{Floratos:1981hs,
+    author = "Floratos, E. G. and Kounnas, C. and Lacaze, R.",
+    title = "{Higher Order QCD Effects in Inclusive Annihilation and Deep Inelastic Scattering}",
+    reportNumber = "LPTENS-81-3",
+    doi = "10.1016/0550-3213(81)90434-X",
+    journal = "Nucl. Phys. B",
+    volume = "192",
+    pages = "417--462",
+    year = "1981"
+}
+
+@article{Moch:2015usa,
+    author = "Moch, S. and Vermaseren, J. A. M. and Vogt, A.",
+    title = "{On \ensuremath{\gamma}5 in higher-order QCD calculations and the NNLO evolution of the polarized valence distribution}",
+    eprint = "1506.04517",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-ph",
+    reportNumber = "DESY-15-061, NIKHEF-15-018, LTH-1042",
+    doi = "10.1016/j.physletb.2015.07.027",
+    journal = "Phys. Lett. B",
+    volume = "748",
+    pages = "432--438",
+    year = "2015"
+}
+
+@article{Moch:2014sna,
+    author = "Moch, S. and Vermaseren, J. A. M. and Vogt, A.",
+    title = "{The Three-Loop Splitting Functions in QCD: The Helicity-Dependent Case}",
+    eprint = "1409.5131",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-ph",
+    reportNumber = "DESY-14-157, NIKHEF-14-033, LTH-1023",
+    doi = "10.1016/j.nuclphysb.2014.10.016",
+    journal = "Nucl. Phys. B",
+    volume = "889",
+    pages = "351--400",
+    year = "2014"
+}

doc/source/theory/N3LO_ad.rst

@@ -4,7 +4,7 @@ N3LO Anomalous Dimensions
 The |N3LO| |QCD| anomalous dimensions :math:`\gamma^{(3)}` are not yet fully known,
 since they rely on the calculation of 4-loop |DIS| integrals.
 Moreover, the analytical structure of these function is already known to be complicated
-since in Mellin space it will included harmonics sum up to weight 7, for which an
+since in Mellin space it include harmonics sum up to weight 7, for which an
 analytical expression is not available.
 
 Here, we describe the various assumptions and limits used in order to reconstruct a parametrization
@@ -331,21 +331,28 @@ final reduced sets of candidates.
         *   - :math:`f_4(N)`
             - :math:`\frac{1}{N-1},\ \frac{1}{N^4},\ \frac{1}{N^3},\ \frac{1}{N^2},\ \frac{1}{N},\ \frac{1}{(N+1)^3},\ \frac{1}{(N+1)^2},\ \frac{1}{N+1},\ \frac{1}{N+2},\ \frac{S_1(N-2)}{N},\ \mathcal{M}[\ln^3(1-x)],\ \mathcal{M}[\ln^2(1-x)], \frac{S_1(N)}{N},\ \frac{S_1^2(N)}{N}`
 
-    Note that this table refers only to the :math:`n_f^0` part where we have assumed no violation of the scaling with :math:`\gamma_{gg}`
+    Note that this table refers only to the :math:`n_f^0` part where we assume no violation of the scaling with :math:`\gamma_{gg}`
     also for the |NLL| term, to help the convergence. We expect that any possible deviation can be parametrized as a shift in he |NNLL| terms
     and in the |NLL| :math:`n_f^1` which are free to vary independently.
+    Furthermore for the part :math:`\propto n_f^2` we adopt a slightly different
+    basis to account fot the fact that the leading
+    contribution for the pole at :math:`N=1` is :math:`\frac{1}{(N-1)^2}`.
 
     .. list-table::  :math:`\gamma_{qg}^{(3)}` parametrization basis
         :align: center
 
         *   - :math:`f_1(N)`
-            - :math:`\frac{1}{(N-1)^2}`
+            - :math:`\frac{1}{(N-1)}`
         *   - :math:`f_2(N)`
             - :math:`\frac{S_1^3(N)}{N}`
         *   - :math:`f_3(N)`
-            - :math:`\frac{1}{N^4},\ \mathcal{M}[\ln^3(1-x)]`
+            - :math:`\frac{1}{N^3},\ \mathcal{M}[\ln^3(1-x)]`
         *   - :math:`f_4(N)`
-            - :math:`\frac{1}{N-1},\ \frac{1}{N^4},\ \frac{1}{N^3},\ \frac{1}{N^2},\ \frac{1}{N},\ \frac{1}{(N+1)^3},\ \frac{1}{(N+1)^2},\ \frac{1}{N+1},\ \frac{1}{N+2},\ \frac{S_1^2(N)}{N},\ \frac{S_1(N-2)}{N},\ \frac{S_1(N)}{N},\ \mathcal{M}[\ln^2(1-x)],\ \mathcal{M}[\ln^3(1-x)]`
+            - :math:`\frac{1}{N^4},\ \frac{1}{N^3},\ \frac{1}{N^2},\ \frac{1}{N},\ \frac{1}{(N+1)^3},\ \frac{1}{(N+1)^2},\ \frac{1}{N+1},\ \frac{1}{N+2},\ \frac{S_1^2(N)}{N},\ \frac{S_1(N-2)}{N},\ \frac{S_1(N)}{N},\ \mathcal{M}[\ln^2(1-x)],\ \mathcal{M}[\ln^3(1-x)]`
+
+    Also in this case we assume no violation of the scaling with :math:`\gamma_{qq}`
+    for the |NLL| term, which is extracted from using 10 moments :cite:`Falcioni:2023luc`
+    thus is known more precisely. We expect that any possible deviation can be parametrized as a shift in he |NNLL| terms.
 
 Slightly different choices are performed for :math:`\gamma_{qq,ps}^{(3)}` where 10 moments
 are known. In this case we can select a larger number of functions in group 3
@@ -366,6 +373,3 @@ and following :cite:`Falcioni:2023luc` we use:
 
 Note that for :math:`\gamma_{qq,ps},\gamma_{qg}` the parts proportional
 to :math:`n_f^0` are not present.
-Furthermore for the part :math:`\propto n_f^2` in :math:`\gamma_{gq}^{(3)}`
-we adopt a slightly different basis to account fot the fact that the leading
-contribution for the pole at :math:`N=1` is :math:`\frac{1}{(N-1)^2}`.

extras/n3lo_bench/plot_mhou.py

@@ -13,17 +13,18 @@
 plt.style.use(utils.load_style())
 
 
-def build_gamma(g, order, q2, fact_sale, nf):
-    a_s = compute_a_s(q2, fact_scale=fact_sale, order=(order + 1, 0), nf=nf)
+def build_gamma(g, order, q2, xif2, nf):
+    a_s = compute_a_s(q2, xif2=xif2, order=(order + 1, 0), nf=nf)
     gamma = 0
     for pto in range(order + 1):
         gamma += g[:, pto] * a_s ** (pto + 1)
     return gamma
 
 
 def compute_mhou(g_min, g_max, g_central, q2, xif2_low, xif2_hig, order, nf):
-    gamma_low = build_gamma(g_min, order, q2 / xif2_low, q2, nf)
-    gamma_hig = build_gamma(g_max, order, q2 / xif2_hig, q2, nf)
+    # TODO: here you need two different instances of strong coupling
+    gamma_low = build_gamma(g_min, order, q2, xif2_low, nf)
+    gamma_hig = build_gamma(g_max, order, q2, xif2_hig, nf)
     delta_low = gamma_low - g_central
     delta_hig = gamma_hig - g_central
     mhou = 0.5 * np.sqrt(delta_low**2 + delta_hig**2)
@@ -43,11 +44,11 @@ def plot_ad(entry, q2=None, nf=4, logscale=True, plot_totu=True, plot_scaling=Fa
     g_low = splitting_function(entry, grid, nf, L=np.log(xif2_low))
     g_hig = splitting_function(entry, grid, nf, L=np.log(xif2_hig))
 
-    a_s_n3lo = compute_a_s(q2, fact_scale=q2, order=(4, 0), nf=nf)
-    g_lo = build_gamma(g, 0, q2, q2, nf)
-    g_nlo = build_gamma(g, 1, q2, q2, nf)
-    g_nnlo = build_gamma(g, 2, q2, q2, nf)
-    g_n3lo = build_gamma(g, 3, q2, q2, nf)
+    a_s_n3lo = compute_a_s(q2, xif2=1, order=(4, 0), nf=nf)
+    g_lo = build_gamma(g, 0, q2, 1, nf)
+    g_nlo = build_gamma(g, 1, q2, 1, nf)
+    g_nnlo = build_gamma(g, 2, q2, 1, nf)
+    g_n3lo = build_gamma(g, 3, q2, 1, nf)
 
     g_n3lo_var = []
     if plot_totu:
@@ -181,10 +182,10 @@ def plot_ad_ratio(entry, q2=None, nf=4, plot_totu=False):
     g_low = splitting_function(entry, grid, nf, L=np.log(xif2_low))
     g_hig = splitting_function(entry, grid, nf, L=np.log(xif2_hig))
 
-    a_s_n3lo = compute_a_s(q2, fact_scale=q2, order=(4, 0), nf=nf)
-    g_nlo = build_gamma(g, 1, q2, q2, nf)
-    g_nnlo = build_gamma(g, 2, q2, q2, nf)
-    g_n3lo = build_gamma(g, 3, q2, q2, nf)
+    a_s_n3lo = compute_a_s(q2, xif2=1, order=(4, 0), nf=nf)
+    g_nlo = build_gamma(g, 1, q2, 1, nf)
+    g_nnlo = build_gamma(g, 2, q2, 1, nf)
+    g_n3lo = build_gamma(g, 3, q2, 1, nf)
 
     g_n3lo_var = []
     if plot_totu:
@@ -282,4 +283,4 @@ def plot_ad_ratio(entry, q2=None, nf=4, plot_totu=False):
 
         # log plots
         x_grid = lambertgrid(80, x_min=1e-7)
-        # plot_ad(k, q2=q2, nf=nf, plot_scaling=True, plot_totu=plot_totu)
+        plot_ad(k, q2=q2, nf=nf, plot_scaling=True, plot_totu=plot_totu)

extras/n3lo_bench/plot_msht.py

@@ -14,10 +14,10 @@
 
 
 n3lo_vars_dict = {
-    "qq": 17,
+    "qq": 8,
     "gg": 18,
     "gq": 24,
-    "qg": 21,
+    "qg": 20,
 }
 
 
@@ -201,4 +201,4 @@ def plot_ad(
 
         # log plots
         x_grid = lambertgrid(60, x_min=1e-7)
-        plot_ad(k, plot_scaling=True, q2=38.5, nf=5)
+        plot_ad(k, plot_scaling=True)

extras/n3lo_bench/splitting_function_utils.py

@@ -98,11 +98,11 @@ def splitting_function(
     return np.array(gamma_x).T
 
 
-def compute_a_s(q2=None, fact_scale=None, nf=None, order=(4, 0)):
+def compute_a_s(q2=None, xif2=1.0, nf=None, order=(4, 0)):
     if q2 is not None:
-        fact_scale = q2 if fact_scale is None else fact_scale
+        fact_scale = q2 * xif2
         ref = types.CouplingsRef(
-            alphas=types.FloatRef(value=0.1181, scale=91.00),
+            alphas=types.FloatRef(value=0.1181, scale=91.00 * np.sqrt(xif2)),
             alphaem=types.FloatRef(value=0.007496, scale=np.nan),
             max_num_flavs=6,
             num_flavs_ref=5,