N3LO msbar and alpha_s

benchmarks/apfel_bench.py

@@ -94,8 +94,9 @@ def benchmark_msbar(self, pto):
         MSbar heavy quark mass scheme
         when  passing kthr != 1 both apfel and eko use ``kThr * msbar``,
         as thr scale, where ``msbar`` is the usual ms_bar solution.
-        However apfel and eko mange the alpha_s thr differently
-        (apfel uses the given mass parameters as thr), so the
+        However apfel and eko manage the alpha_s thr differently:
+        apfel uses the given mass parameters as thr multiplied by the kthr,
+        while in eko only the already computed thr matters, so the
         benchmark is not a proper comparison with this option allowed.
         """
         th = self.vfns_theory.copy()
@@ -179,4 +180,4 @@ def benchmark_sv(self, pto):
     # obj.benchmark_plain(1)
     # obj.benchmark_sv(1)
     # obj.benchmark_kthr(2)
-    obj.benchmark_msbar(0)
+    obj.benchmark_msbar(2)

doc/source/refs.bib

@@ -419,3 +419,30 @@ @article{Bertone:2013vaa
     pages         = {1647--1668},
     year          = {2014}
 }
+
+@article{Vermaseren:1997fq,
+    author = "Vermaseren, J. A. M. and Larin, S. A. and van Ritbergen, T.",
+    title = "{The four loop quark mass anomalous dimension and the invariant quark mass}",
+    eprint = "hep-ph/9703284",
+    archivePrefix = "arXiv",
+    reportNumber = "UM-TH-97-03, NIKHEF-97-012",
+    doi = "10.1016/S0370-2693(97)00660-6",
+    journal = "Phys. Lett. B",
+    volume = "405",
+    pages = "327--333",
+    year = "1997"
+}
+
+@article{Liu:2015fxa,
+    author = "Liu, Tao and Steinhauser, Matthias",
+    title = "{Decoupling of heavy quarks at four loops and effective Higgs-fermion coupling}",
+    eprint = "1502.04719",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-ph",
+    reportNumber = "TTP15-05",
+    doi = "10.1016/j.physletb.2015.05.023",
+    journal = "Phys. Lett. B",
+    volume = "746",
+    pages = "330--334",
+    year = "2015"
+}

doc/source/theory/pQCD.rst

@@ -23,13 +23,29 @@ We implement two different strategies to solve the |RGE|:
 - ``method="expanded"``: using approximate solutions:
 
 .. math ::
-    a^{\text{LO}}_s(\mu_R^2) &= \frac{a_s(\mu_0^2)}{1 + a_s(\mu_0^2) \beta_0 \ln(\mu_R^2/\mu_0^2)} \\
-    a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2) &= a^{\text{LO}}_s(\mu_R^2)-b_1 \left[a^{\text{LO}}_s(\mu_R^2)\right]^2 \ln\left(1+a_s(\mu_0^2) \beta_0 \ln(\mu_R^2/\mu_0^2)\right) \\
+    a^{\text{LO}}_s(\mu_R^2) &= \frac{a_s(\mu_0^2)}{1 + a_s(\mu_0^2) \beta_0 L_{\mu}} \\
+    a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2) &= a^{\text{LO}}_s(\mu_R^2)-b_1 \left[a^{\text{LO}}_s(\mu_R^2)\right]^2 \ln\left(1+a_s(\mu_0^2) \beta_0 L_{\mu}\right) \\
     a^{\text{NNLO}}_{s,\text{exp}}(\mu_R^2) &= a^{\text{LO}}_s(\mu_R^2)\left[1 + a^{\text{LO}}_s(\mu_R^2)\left(a^{\text{LO}}_s(\mu_R^2) - a_s(\mu_0^2)\right)(b_2 - b_1^2) \right.\\
-                                        & \hspace{60pt} \left. + a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2) b_1 \ln\left(a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2)/a_s(\mu_0^2)\right)\right]
+                                        & \hspace{60pt} \left. + a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2) b_1 \ln\left(a^{\text{NLO}}_{s,\text{exp}}(\mu_R^2)/a_s(\mu_0^2)\right)\right] \\
+    a^{\text{N3LO}}_{s,\text{exp}}(\mu_R^2) &= a^{\text{NNLO}}_s(\mu_R^2) + \frac{a^{\text{LO}}_s(\mu_R^2)^4}{2 b_0^3} \left\{ \right. \\
+                & -2 b_1^3 L_{0}^3 + 5 b_1^3 L_{\text{LO}}^2 + 2 b_1^3  L_{\text{LO}}^3 + b_1^3 L_{0}^2 \left(5 + 6  L_{\text{LO}} \right) \\
+                & + 2 b_0 b_1  L_{\text{LO}} \left[ b_2 + 2 \left(b_1^2 - b_0 b_2 \right) L_{\mu} a_s(\mu_0^2) \right] \\
+                & - b_0^2 L_{\mu} a_s(\mu_0^2) \left[ -2 b_1 b_2 + 2 b_0 b_3 + \left( b_1^3 - 2 b_0 b_1 b_2 + b_0^2 b_3 \right) L_{\mu} a_s(\mu_0^2) \right] \\
+                & - 2 b_1 L_{0} \left[ 5 b_1^2  L_{\text{LO}} + 3 b_1^2  L_{\text{LO}}^2 + b_0 \left[b_2 + 2 \left(b_1^2 - b_0 b_2\right) L_{\mu} a_s(\mu_0^2)\right] \right ] \\
+                & \left. \right\}
+
+being:
+
+.. math ::
+    L_{\mu} &= \ln(\mu_R^2/\mu_0^2) \\
+    L_{0} &= \ln(a_s(\mu_0^2)) \\
+    L_{\text{LO}} &= \ln(a^{\text{LO}}_s(\mu_R^2)) \\
 
 When the renormalization scale crosses a flavor threshold matching conditions
 have to be applied :cite:`Schroder:2005hy,Chetyrkin:2005ia`.
+In particular, the matching involved in the change from :math:`n_f` to :math:`n_f-1` schemes
+is presented in equation 3.1 of :cite:`Schroder:2005hy` for |MSbar| masses, while the
+same expression for POLE masses is reported in Appendix A.
 
 Splitting Functions
 -------------------
@@ -97,7 +113,7 @@ For each heavy quark :math:`h` we solve for :math:`m_h`:
 where the evolved |MSbar| mass is given by:
 
 .. math ::
-    m_{\overline{MS},h}(\mu^2) = m_{h,0} \int_{a_s(\mu_{h,0}^2)}^{a_s(\mu^2)} \frac{\gamma(a_s)}{\beta(a_s)} d a_s
+    m_{\overline{MS},h}(\mu^2) = m_{h,0} \exp \left[ - \int_{a_s(\mu_{h,0}^2)}^{a_s(\mu^2)} \frac{\gamma_m(a_s)}{\beta(a_s)} d a_s \right ]
 
 and :math:`m_{h,0}` is the given initial condition at the scale
 :math:`\mu_{h,0}`. Here there is a subtle complication since the solution
@@ -127,24 +143,28 @@ In doing so EKO takes advantage of the monotony of the |RGE| solution
 
 Now, being able to evaluate :math:`a_s(\mu_{h,0}^2)`, there are two ways of
 solving the previous integral and finally compute the evolved
-:math:`m_{\overline{MS},h}`. In fact, the function :math:`\gamma(a_s)` is the
+:math:`m_{\overline{MS},h}`. In fact, the function :math:`\gamma_m(a_s)` is the
 anomalous QCD mass dimension and, as the :math:`\beta` function, it can be evaluated
-perturbatively in :math:`a_s` up to :math:`\mathcal{O}(a_s^3)`:
+perturbatively in :math:`a_s` up to :math:`\mathcal{O}(a_s^4)`:
 
 .. math ::
-    \gamma(a_s) &= - \sum\limits_{n=0} \gamma_n a_s^{n+1} \\
+    \gamma_m(a_s) &= \sum\limits_{n=0} \gamma_{m,n} a_s^{n+1} \\
 
-Even here it is useful to define :math:`c_k = \gamma_k/\beta_0, k>0`.
+Even here it is useful to define :math:`c_k = \gamma_{m,k}/\beta_0, k \ge 0`.
 
 Therefore the two solution strategies are:
 
 - ``method = "exact"``: the integral is solved exactly using the expression of
-  :math:`\beta,\gamma` up to the specified perturbative order
+  :math:`\beta,\gamma_m` up to the specified perturbative order
 - ``method = "expanded"``: the integral is approximate by the following expansion:
 
 .. math ::
     m_{\overline{MS},h}(\mu^2) & = m_{h,0} \left ( \frac{a_s(\mu^2)}{a_s(\mu_{h,0}^2)} \right )^{c_0} \frac{j_{exp}(a_s(\mu^2))}{j_{exp}(a_s(\mu_{h,0}^2))} \\
-    j_{exp}(a_s) &= 1 + a_s \left [ c_1 - b_1 c_0 \right ] + \frac{a_s^2}{2} \left [c_2 - c_1 b_1 - b_2 c_0 + b_1^2 c_0 + (c_1 - b_1 c_0)^2 \right]
+    j_{exp}(a_s) &= 1 + a_s \left [ c_1 - b_1 c_0 \right ] \\
+                 & + \frac{a_s^2}{2} \left [c_2 - c_1 b_1 - b_2 c_0 + b_1^2 c_0 + (c_1 - b_1 c_0)^2 \right] \\
+                 & + \frac{a_s^3}{6} [ -2 b_3 c_0 - b_1^3 c_0 (1 + c_0) (2 + c_0) - 2 b_2 c_1 \\
+                 & - 3 b_2 c_0 c_1 + b_1^2 (2 + 3 c_0 (2 + c_0)) c_1 + c_1^3 + 3 c_1 c_2 \\
+                 & + b_1 (b_2 c_0 (4 + 3 c_0) - 3 (1 + c_0) c_1^2 - (2 + 3 c_0) c_2) + 2 c_3 ]
 
 
 The procedure is iterated on all the heavy quarks, updating the temporary instance
@@ -154,10 +174,11 @@ To find coherent solutions and perform the mass running in the correct patches i
 is necessary to always start computing the mass scales closer to :math:`\mu_{ref}`.
 
 Eventually, to ensure that the threshold values are properly set, we add a
-consistency check, asserting:
+consistency check, asserting that the :math:`m_{\overline{MS},h}` are properly sorted.
 
-.. math ::
-    m_{\overline{MS},h} (m_h) \leq m_{\overline{MS},h+1} (m_h)
+Note that also for |MSbar| mass running when the heavy threshold scales are
+crossed we need to apply non trivial matching from order
+:math:`\mathcal{O}(a_s^2)` as described here :cite:`Liu:2015fxa`.
 
 We provide the following as an illustrative example of how this procedure works:
 when the strong coupling is given with boundary condition :math:`\alpha_s(\mu_{ref}=91, n_{f_{ref}}=5)`

src/eko/anomalous_dimensions/harmonics.py

@@ -13,10 +13,11 @@
 zeta2 = scipy.special.zeta(2)
 zeta3 = scipy.special.zeta(3)
 zeta4 = scipy.special.zeta(4)
+zeta5 = scipy.special.zeta(5)
 
 
 @nb.njit("c16(c16,u1)", cache=True)
-def cern_polygamma(Z, K: int):  # pylint: disable=all
+def cern_polygamma(Z, K):  # pylint: disable=all
     """
     Computes the polygamma functions :math:`\\psi_k(z)`.
 

src/eko/anomalous_dimensions/lo.py

@@ -32,7 +32,7 @@ def gamma_ns_0(N, s1):
 
 
 @nb.njit("c16(c16,u1)", cache=True)
-def gamma_qg_0(N, nf: int):
+def gamma_qg_0(N, nf):
     """
     Computes the leading-order quark-gluon anomalous dimension
 
@@ -78,7 +78,7 @@ def gamma_gq_0(N):
 
 
 @nb.njit("c16(c16,c16,u1)", cache=True)
-def gamma_gg_0(N, s1, nf: int):
+def gamma_gg_0(N, s1, nf):
     """
     Computes the leading-order gluon-gluon anomalous dimension
 
@@ -104,7 +104,7 @@ def gamma_gg_0(N, s1, nf: int):
 
 
 @nb.njit("c16[:,:](c16,c16,u1)", cache=True)
-def gamma_singlet_0(N, s1, nf: int):
+def gamma_singlet_0(N, s1, nf):
     r"""
       Computes the leading-order singlet anomalous dimension matrix
 

src/eko/anomalous_dimensions/nlo.py

@@ -15,7 +15,7 @@
 
 
 @nb.njit("c16(c16,u1)", cache=True)
-def gamma_nsm_1(n, nf: int):
+def gamma_nsm_1(n, nf):
     """
     Computes the |NLO| valence-like non-singlet anomalous dimension.
 
@@ -58,7 +58,7 @@ def gamma_nsm_1(n, nf: int):
 
 
 @nb.njit("c16(c16,u1)", cache=True)
-def gamma_nsp_1(n, nf: int):
+def gamma_nsp_1(n, nf):
     """
     Computes the |NLO| singlet-like non-singlet anomalous dimension.
 
@@ -99,7 +99,7 @@ def gamma_nsp_1(n, nf: int):
 
 
 @nb.njit("c16(c16,u1)", cache=True)
-def gamma_ps_1(n, nf: int):
+def gamma_ps_1(n, nf):
     """
     Computes the |NLO| pure-singlet quark-quark anomalous dimension.
 
@@ -126,7 +126,7 @@ def gamma_ps_1(n, nf: int):
 
 
 @nb.njit("c16(c16,u1)", cache=True)
-def gamma_qg_1(n, nf: int):
+def gamma_qg_1(n, nf):
     """
     Computes the |NLO| quark-gluon singlet anomalous dimension.
 
@@ -159,7 +159,7 @@ def gamma_qg_1(n, nf: int):
 
 
 @nb.njit("c16(c16,u1)", cache=True)
-def gamma_gq_1(n, nf: int):
+def gamma_gq_1(n, nf):
     """
     Computes the |NLO| gluon-quark singlet anomalous dimension.
 
@@ -195,7 +195,7 @@ def gamma_gq_1(n, nf: int):
 
 
 @nb.njit("c16(c16,u1)", cache=True)
-def gamma_gg_1(n, nf: int):
+def gamma_gg_1(n, nf):
     """
     Computes the |NLO| gluon-gluon singlet anomalous dimension.
 
@@ -233,7 +233,7 @@ def gamma_gg_1(n, nf: int):
 
 
 @nb.njit("c16[:,:](c16,u1)", cache=True)
-def gamma_singlet_1(N, nf: int):
+def gamma_singlet_1(N, nf):
     r"""
       Computes the next-leading-order singlet anomalous dimension matrix
 

src/eko/anomalous_dimensions/nnlo.py

@@ -17,7 +17,7 @@
 
 
 @nb.njit("c16(c16,u1,c16[:])", cache=True)
-def gamma_nsm_2(n, nf: int, sx):
+def gamma_nsm_2(n, nf, sx):
     """
     Computes the |NNLO| valence-like non-singlet anomalous dimension.
 
@@ -94,7 +94,7 @@ def gamma_nsm_2(n, nf: int, sx):
 
 
 @nb.njit("c16(c16,u1,c16[:])", cache=True)
-def gamma_nsp_2(n, nf: int, sx):
+def gamma_nsp_2(n, nf, sx):
     """
     Computes the |NNLO| singlet-like non-singlet anomalous dimension.
 
@@ -171,7 +171,7 @@ def gamma_nsp_2(n, nf: int, sx):
 
 
 @nb.njit("c16(c16,u1,c16[:])", cache=True)
-def gamma_nsv_2(n, nf: int, sx):
+def gamma_nsv_2(n, nf, sx):
     """
     Computes the |NNLO| valence non-singlet anomalous dimension.
 
@@ -226,7 +226,7 @@ def gamma_nsv_2(n, nf: int, sx):
 
 
 @nb.njit("c16(c16,u1,c16[:])", cache=True)
-def gamma_ps_2(n, nf: int, sx):
+def gamma_ps_2(n, nf, sx):
     """
     Computes the |NNLO| pure-singlet quark-quark anomalous dimension.
 
@@ -298,7 +298,7 @@ def gamma_ps_2(n, nf: int, sx):
 
 
 @nb.njit("c16(c16,u1,c16[:])", cache=True)
-def gamma_qg_2(n, nf: int, sx):
+def gamma_qg_2(n, nf, sx):
     """
     Computes the |NNLO| quark-gluon singlet anomalous dimension.
 
@@ -372,7 +372,7 @@ def gamma_qg_2(n, nf: int, sx):
 
 
 @nb.njit("c16(c16,u1,c16[:])", cache=True)
-def gamma_gq_2(n, nf: int, sx):
+def gamma_gq_2(n, nf, sx):
     """
     Computes the |NNLO| gluon-quark singlet anomalous dimension.
 
@@ -462,7 +462,7 @@ def gamma_gq_2(n, nf: int, sx):
 
 
 @nb.njit("c16(c16,u1,c16[:])", cache=True)
-def gamma_gg_2(n, nf: int, sx):
+def gamma_gg_2(n, nf, sx):
     """
     Computes the |NNLO| gluon-gluon singlet anomalous dimension.
 
@@ -550,7 +550,7 @@ def gamma_gg_2(n, nf: int, sx):
 
 
 @nb.njit("c16[:,:](c16,u1,c16[:])", cache=True)
-def gamma_singlet_2(N, nf: int, sx):
+def gamma_singlet_2(N, nf, sx):
     r"""
       Computes the |NNLO| singlet anomalous dimension matrix
 

src/eko/beta.py

@@ -8,10 +8,11 @@
 import numba as nb
 
 from eko import constants
+from eko.anomalous_dimensions.harmonics import zeta3
 
 
 @nb.njit("f8(u1)", cache=True)
-def beta_0(nf: int):
+def beta_0(nf):
     """
     Computes the first coefficient of the QCD beta function.
 
@@ -32,7 +33,7 @@ def beta_0(nf: int):
 
 
 @nb.njit("f8(u1)", cache=True)
-def beta_1(nf: int):
+def beta_1(nf):
     """
     Computes the second coefficient of the QCD beta function.
 
@@ -57,7 +58,7 @@ def beta_1(nf: int):
 
 
 @nb.njit("f8(u1)", cache=True)
-def beta_2(nf: int):
+def beta_2(nf):
     """
     Computes the third coefficient of the QCD beta function
 
@@ -85,6 +86,33 @@ def beta_2(nf: int):
     return beta_2
 
 
+@nb.njit("f8(u1)", cache=True)
+def beta_3(nf):
+    """
+    Computes the fourth coefficient of the QCD beta function
+
+    Implements Eq. (3.6) of :cite:`Herzog:2017ohr`.
+
+    Parameters
+    ----------
+        nf : int
+            number of active flavors
+
+    Returns
+    -------
+        beta_3 : float
+            fourth coefficient of the QCD beta function :math:`\\beta_3^{n_f}`
+    """
+    beta_3 = (
+        149753.0 / 6.0
+        + 3564.0 * zeta3
+        + nf * (-1078361.0 / 162.0 - 6508.0 / 27.0 * zeta3)
+        + nf**2 * (50065.0 / 162.0 + 6472.0 / 81.0 * zeta3)
+        + 1093.0 / 729.0 * nf**3
+    )
+    return beta_3
+
+
 @nb.njit("f8(u1,u1)", cache=True)
 def beta(k, nf):
     """
@@ -109,8 +137,10 @@ def beta(k, nf):
         beta_ = beta_1(nf)
     elif k == 2:
         beta_ = beta_2(nf)
+    elif k == 3:
+        beta_ = beta_3(nf)
     else:
-        raise ValueError("Beta coefficients beyond NNLO are not implemented!")
+        raise ValueError("Beta coefficients beyond N3LO are not implemented!")
     return beta_