diff --git a/src/eko/evolution_operator/__init__.py b/src/eko/evolution_operator/__init__.py index 323b117b4..5818c90fc 100644 --- a/src/eko/evolution_operator/__init__.py +++ b/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( diff --git a/src/eko/matching_conditions/operator_matrix_element.py b/src/eko/matching_conditions/operator_matrix_element.py index a8ad60109..ed0299294 100644 --- a/src/eko/matching_conditions/operator_matrix_element.py +++ b/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)""" diff --git a/src/eko/scale_variations/expanded.py b/src/eko/scale_variations/expanded.py index 7eb35750e..3ceb5c6d4 100644 --- a/src/eko/scale_variations/expanded.py +++ b/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 diff --git a/tests/eko/evolution_operator/test_init.py b/tests/eko/evolution_operator/test_init.py index dc8ddbd8e..d40fbeb39 100644 --- a/tests/eko/evolution_operator/test_init.py +++ b/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, diff --git a/tests/eko/scale_variations/test_expanded.py b/tests/eko/scale_variations/test_expanded.py index bce253f9c..c1ecc33fd 100644 --- a/tests/eko/scale_variations/test_expanded.py +++ b/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", )