Fix scale variations expanded

src/eko/evolution_operator/__init__.py

@@ -123,7 +123,6 @@ def quad_ker(
     areas,
     as1,
     as0,
-    as_raw,
     nf,
     L,
     ev_op_iterations,
@@ -155,8 +154,6 @@ def quad_ker(
         target coupling value
     as0 : float
         initial coupling value
-    as_raw : float
-        coupling value at the process scale
     nf : int
         number of active flavors
     L : float
@@ -201,7 +198,7 @@ def quad_ker(
         # scale var expanded is applied on the kernel
         if sv_mode == sv.Modes.expanded and not is_threshold:
             ker = np.ascontiguousarray(
-                sv.expanded.singlet_variation(gamma_singlet, as_raw, order, nf, L)
+                sv.expanded.singlet_variation(gamma_singlet, as1, order, nf, L)
             ) @ np.ascontiguousarray(ker)
         ker = select_singlet_element(ker, mode0, mode1)
     else:
@@ -218,9 +215,7 @@ def quad_ker(
             ev_op_iterations,
         )
         if sv_mode == sv.Modes.expanded and not is_threshold:
-            ker = (
-                sv.expanded.non_singlet_variation(gamma_ns, as_raw, order, nf, L) * ker
-            )
+            ker = sv.expanded.non_singlet_variation(gamma_ns, as1, order, nf, L) * ker
 
     # recombine everything
     return np.real(ker * integrand)
@@ -291,7 +286,7 @@ def grid_size(self):
         """Return the grid size."""
         return self.int_disp.xgrid.size
 
-    def mur2_shift(self, q2):
+    def sv_exponentiated_shift(self, q2):
         """Compute shifted renormalization scale.
 
         Parameters
@@ -313,11 +308,16 @@ def a_s(self):
         """Return the computed values for :math:`a_s`."""
         sc = self.managers["strong_coupling"]
         a0 = sc.a_s(
-            self.mur2_shift(self.q2_from), fact_scale=self.q2_from, nf_to=self.nf
+            self.sv_exponentiated_shift(self.q2_from),
+            fact_scale=self.q2_from,
+            nf_to=self.nf,
+        )
+        a1 = sc.a_s(
+            self.sv_exponentiated_shift(self.q2_to),
+            fact_scale=self.q2_to,
+            nf_to=self.nf,
         )
-        a1 = sc.a_s(self.mur2_shift(self.q2_to), fact_scale=self.q2_to, nf_to=self.nf)
-        a_raw = sc.a_s(self.q2_to, fact_scale=self.q2_to, nf_to=self.nf)
-        return (a0, a1, a_raw)
+        return (a0, a1)
 
     @property
     def labels(self):
@@ -375,7 +375,6 @@ def quad_ker(self, label, logx, areas):
             areas=areas,
             as1=self.a_s[1],
             as0=self.a_s[0],
-            as_raw=self.a_s[2],
             nf=self.nf,
             L=np.log(self.xif2),
             ev_op_iterations=self.config["ev_op_iterations"],
@@ -473,8 +472,8 @@ def compute(self):
         logger.info(
             "%s: µ_R^2 distance: %e -> %e",
             self.log_label,
-            self.mur2_shift(self.q2_from),
-            self.mur2_shift(self.q2_to),
+            self.sv_exponentiated_shift(self.q2_from),
+            self.sv_exponentiated_shift(self.q2_to),
         )
         if self.sv_mode != sv.Modes.unvaried:
             logger.info(

src/eko/matching_conditions/operator_matrix_element.py

@@ -155,7 +155,7 @@ def build_ome(A, matching_order, a_s, backward_method):
 def quad_ker(
     u, order, mode0, mode1, is_log, logx, areas, a_s, nf, L, backward_method, is_msbar
 ):
-    """
+    r"""
     Raw kernel inside quad
 
     Parameters
@@ -366,7 +366,9 @@ def a_s(self):
         Note that here you need to use :math:`a_s^{n_f+1}`
         """
         sc = self.managers["strong_coupling"]
-        return sc.a_s(self.mur2_shift(self.q2_from), self.q2_from, nf_to=self.nf + 1)
+        return sc.a_s(
+            self.sv_exponentiated_shift(self.q2_from), self.q2_from, nf_to=self.nf + 1
+        )
 
     def compute(self):
         """Compute the actual operators (i.e. run the integrations)"""

src/eko/scale_variations/expanded.py

@@ -1,4 +1,4 @@
-r"""This module contains the scale variation operator for the expanded scheme (``ModSV=expanded``).
+r"""Scale variation operator for the expanded scheme (``ModSV=expanded``).
 
 The expressions can be obtained using Eqs. (3.33) and (3.38) of :cite:`AbdulKhalek:2019ihb`.
 Be aware that corresponding the signs of the ingredients there are a number of differences.
@@ -14,7 +14,7 @@
 
 @nb.njit(cache=True)
 def variation_as1(gamma, L):
-    r"""Computes the |NLO| anomalous dimension variation.
+    r"""Compute the |NLO| anomalous dimension variation.
 
     Parameters
     ----------
@@ -33,7 +33,11 @@ def variation_as1(gamma, L):
 
 @nb.njit(cache=True)
 def variation_as2(gamma, L, beta0, g0e2):
-    r"""Computes the |NNLO| anomalous dimension variation.
+    r"""Compute the |NNLO| anomalous dimension variation.
+
+    These kernels are meant to be used with alpha_s evaluated at the
+    factorization scale. If one expresses everything in terms of
+    alpha_s evaluated at the process scale, the sign of beta*gamma0 flips.
 
     Parameters
     ----------
@@ -56,7 +60,7 @@ def variation_as2(gamma, L, beta0, g0e2):
 
 @nb.njit(cache=True)
 def variation_as3(gamma, L, beta0, beta1, g0e2, g0e3, g1g0, g0g1):
-    r"""Computes the |N3LO| anomalous dimension variation.
+    r"""Compute the |N3LO| anomalous dimension variation.
 
     Parameters
     ----------
@@ -117,14 +121,14 @@ def non_singlet_variation(gamma, a_s, order, nf, L):
     """
     sv_ker = 1.0
     if order[0] >= 2:
-        sv_ker += a_s * variation_as1(gamma, L)
+        sv_ker -= a_s * variation_as1(gamma, L)
     if order[0] >= 3:
         beta0 = beta.beta_qcd_as2(nf)
-        sv_ker += a_s**2 * variation_as2(gamma, L, beta0, gamma[0] ** 2)
+        sv_ker -= a_s**2 * variation_as2(gamma, L, beta0, gamma[0] ** 2)
     if order[0] >= 4:
         beta1 = beta.beta_qcd((3, 0), nf)
         g0g1 = gamma[0] * gamma[1]
-        sv_ker += a_s**3 * variation_as3(
+        sv_ker -= a_s**3 * variation_as3(
             gamma, L, beta0, beta1, gamma[0] ** 2, gamma[0] ** 3, g0g1, g0g1
         )
     return sv_ker
@@ -155,18 +159,18 @@ def singlet_variation(gamma, a_s, order, nf, L):
     sv_ker = np.eye(2, dtype=np.complex_)
     gamma = np.ascontiguousarray(gamma)
     if order[0] >= 2:
-        sv_ker += a_s * variation_as1(gamma, L)
+        sv_ker -= a_s * variation_as1(gamma, L)
     if order[0] >= 3:
         beta0 = beta.beta_qcd_as2(nf)
         gamma0e2 = gamma[0] @ gamma[0]
-        sv_ker += a_s**2 * variation_as2(gamma, L, beta0, gamma0e2)
+        sv_ker -= a_s**2 * variation_as2(gamma, L, beta0, gamma0e2)
     if order[0] >= 4:
         beta1 = beta.beta_qcd((3, 0), nf)
         gamma0e3 = gamma0e2 @ gamma[0]
         # here the product is not commutative
         g1g0 = gamma[1] @ gamma[0]
         g0g1 = gamma[0] @ gamma[1]
-        sv_ker += a_s**3 * variation_as3(
+        sv_ker -= a_s**3 * variation_as3(
             gamma, L, beta0, beta1, gamma0e2, gamma0e3, g1g0, g0g1
         )
     return sv_ker

tests/eko/evolution_operator/test_init.py

@@ -37,7 +37,6 @@ def test_quad_ker(monkeypatch):
             areas=np.zeros(3),
             as1=1,
             as0=2,
-            as_raw=1,
             nf=3,
             L=0,
             ev_op_iterations=0,
@@ -57,7 +56,6 @@ def test_quad_ker(monkeypatch):
             areas=np.zeros(3),
             as1=1,
             as0=2,
-            as_raw=1,
             nf=3,
             L=0,
             ev_op_iterations=0,
@@ -77,7 +75,6 @@ def test_quad_ker(monkeypatch):
             areas=np.zeros(3),
             as1=1,
             as0=2,
-            as_raw=1,
             nf=3,
             L=0,
             ev_op_iterations=0,
@@ -99,7 +96,6 @@ def test_quad_ker(monkeypatch):
                 areas=np.zeros(3),
                 as1=1,
                 as0=2,
-                as_raw=1,
                 nf=3,
                 L=0,
                 ev_op_iterations=0,
@@ -121,7 +117,6 @@ def test_quad_ker(monkeypatch):
         areas=np.zeros(3),
         as1=1,
         as0=2,
-        as_raw=1,
         nf=3,
         L=0,
         ev_op_iterations=0,
@@ -189,7 +184,7 @@ def test_exponentiated(self, theory_ffns, operator_card, tmp_path):
         g = r.op_grid
         # setup objs
         o = Operator(g.config, g.managers, 3, 2.0, 10.0)
-        np.testing.assert_allclose(o.mur2_shift(40.0), 10.0)
+        np.testing.assert_allclose(o.sv_exponentiated_shift(40.0), 10.0)
         o.compute()
         self.check_lo(o)
 
@@ -314,7 +309,7 @@ def quad_ker_pegasus(
     mode1 = 0
     method = ""
     logxs = np.log(int_disp.xgrid.raw)
-    as_raw = a1 = 1
+    a1 = 1
     a0 = 2
     nf = 3
     L = 0
@@ -335,7 +330,6 @@ def quad_ker_pegasus(
                     bf.areas_representation,
                     a1,
                     a0,
-                    as_raw,
                     nf,
                     L,
                     ev_op_iterations,

tests/eko/scale_variations/test_expanded.py

@@ -2,7 +2,7 @@
 
 from eko import basis_rotation as br
 from eko.anomalous_dimensions import gamma_ns, gamma_singlet
-from eko.beta import beta_qcd_as2
+from eko.beta import beta_qcd_as2, beta_qcd_as3
 from eko.kernels import non_singlet, singlet
 from eko.scale_variations import Modes, expanded, exponentiated
 
@@ -43,13 +43,12 @@ def test_singlet_sv_dispacher():
 def test_scale_variation_a_vs_b():
     r"""
     Test ``ModSV=exponentiated`` kernel vs ``ModSV=expanded``.
-    We test that the quantity :math:`(ker_A - ker_B)/ker_{unv}` depends
-    only on the accuracy in the :math:`\alpha_s` expansion
-    and not in the size of the `fact_to_ren` itself.
-    However in our implementation the exponentiated mode depends on
-    the actual value of a0 and a1, since the evolution integral in :math:`\alpha_s`
-    is evaluated. Thus this test ratio :math:`(ker_A - ker_B)/ker_{unv}`
-    still contains a dependency on `fact_to_ren`.
+    We test that the quantity :math:`(ker_A - ker_B)/ker_{unv}`
+
+    Note this is NOT the real difference between scheme expanded
+    and exponentiated since here we don't take into account the
+    shifts in path length and :math:`\alpha_s` values.
+    The real difference is always higher order.
     """
     nf = 5
     n = 10
@@ -64,20 +63,50 @@ def scheme_diff(g, k, pto, is_singlet):
         the accuracy of singlet quantities is slightly worse.
         """
         if pto[0] >= 2:
-            diff = g[0] * k * a0
+            diff = g[0] * k * a0 - 2 * a1 * k * g[0]
         if pto[0] >= 3:
             b0 = beta_qcd_as2(nf)
             g02 = g[0] @ g[0] if is_singlet else g[0] ** 2
-            diff += a0**2 * g[1] * k - k**2 * (
-                1 / 2 * a0**2 * b0 * g[0] + a1 * a0 * g02 - 1 / 2 * a0**2 * g02
+            diff += (
+                -2 * a1**2 * g[1] * k
+                + a0**2 * g[1] * k
+                + a1**2 * b0 * g[0] * k**2
+                - 0.5 * a0**2 * b0 * g[0] * k**2
+                - a1 * a0 * g02 * k**2
+                + 0.5 * a0**2 * g02 * k**2
+                + a1**2 * g02 * k**2
+            )
+        if pto[0] >= 4:
+            b1 = beta_qcd_as3(nf)
+            g0g1 = g[0] @ g[1] if is_singlet else g[0] * g[1]
+            g03 = g02 @ g[0] if is_singlet else g02 * g[0]
+            diff += (
+                a0**3 * g[2] * k
+                - 2 * a1**3 * g[2] * k
+                - 1 / 2 * a0**3 * b1 * g[0] * k**2
+                + a1**3 * b1 * g[0] * k**2
+                - a0**3 * b0 * g[1] * k**2
+                + 2 * a1**3 *b0 * g[1] * k**2
+                + a0**3 * g0g1 * k**2
+                - a0**2 * a1 * g0g1 * k**2
+                - a0 * a1**2 * g0g1 * k**2
+                + 2 * a1**3 * g0g1 * k**2
+                + 1 / 3 * a0**3 * b0**2 * g[0] * k**3
+                - 2 / 3 * a1**3 * b0**2 * g[0] * k**3
+                - 1 / 2 * a0**3 * b0 * g02 * k**3
+                + 1 / 2 * a0**2 * a1 * b0 * g02 * k**3
+                + 1 / 2 * a0 * a1**2 * b0 * g02 * k**3
+                - a1**3 * b0 * g02 * k**3
+                + 1 / 6 * a0**3 * g03 * k**3
+                - 1 / 2 * a0**2 * a1 * g03 * k**3
+                + 1 / 2 * a0 * a1**2 * g03 * k**3
+                - 1 / 3 * a1**3 * g03 * k**3
             )
         return diff
 
-    # TODO: perform this test also at N3LO, once evolution kernels
-    # will be implemented
 
     for L in [np.log(0.5), np.log(2)]:
-        for order in [(2, 0), (3, 0)]:
+        for order in [(2, 0), (3, 0), (4,0)]:
             # Non singlet kernels
             gns = gamma_ns(order, br.non_singlet_pids_map["ns+"], n, nf)
             ker = non_singlet.dispatcher(
@@ -89,7 +118,12 @@ def scheme_diff(g, k, pto, is_singlet):
             )
             ker_b = ker * expanded.non_singlet_variation(gns, a1, order, nf, L)
             ns_diff = scheme_diff(gns, L, order, False)
-            np.testing.assert_allclose((ker_a - ker_b) / ker, ns_diff, atol=1e-3)
+            np.testing.assert_allclose(
+                (ker_a - ker_b) / ker,
+                ns_diff,
+                atol=1e-3,
+                err_msg=f"L={L},order={order},non-singlet",
+            )
 
             # Singlet kernels
             gs = gamma_singlet(order, n, nf)
@@ -103,5 +137,8 @@ def scheme_diff(g, k, pto, is_singlet):
             ker_b = ker @ expanded.singlet_variation(gs, a1, order, nf, L)
             s_diff = scheme_diff(gs, L, order, True)
             np.testing.assert_allclose(
-                (ker_a - ker_b) @ np.linalg.inv(ker), s_diff, atol=5e-3
+                (ker_a - ker_b) @ np.linalg.inv(ker),
+                s_diff,
+                atol=5e-3,
+                err_msg=f"L={L},order={order},singlet",
             )