On the perturbative renormalization of four-quark operators for new physics

Abstract

We discuss the renormalization properties of the full set of $\Delta F=2$ operators involved in BSM processes, including the definition of RGI versions of operators that exhibit mixing under RG transformations. As a first step for a fully non-perturbative determination of the scale-dependent renormalization factors and their runnings, we introduce a family of appropriate Schrödinger Functional schemes, and study them in perturbation theory. This allows, in particular, to determine the NLO anomalous dimensions of all $\Delta F=1,2$ operators in these schemes. Finally, we discuss the systematic uncertainties related to the use of NLO perturbation theory for the RG running of four-quark operators to scales in the GeV range, in both our SF schemes and standard ${\overline{\mathrm{MS}}}$ and RI-MOM schemes. Large truncation effects are found for some of the operators considered.

1 Introduction

Hadronic matrix elements of four-quark operators play an important rôle in the study of flavor physics within the Standard Model (SM), as well as in searches for new physics. In particular, they are essential to the study of CP-violation in the hadron sector in both the SM and beyond-the-SM (BSM) models, where they parametrize the effect of new interactions. A key ingredient of these studies is the renormalization of the operators, including their renormalization group (RG) running from low-energy hadronic scales $\mathcal {O}(\Lambda _\mathrm{\scriptscriptstyle QCD})$ to the high-energy electroweak or new physics scales, where contact with the fundamental underlying theory is made.

In this paper we prepare the ground for a full non-perturbative computation of the low-energy renormalization and RG running of all possible four-quark operators with net flavor change, by introducing appropriate Schrödinger Functional (SF) renormalization schemes. In order to connect them with standard ${\overline{\mathrm{MS}}}$ schemes at high energies, as well as with renormalization group invariant (RGI) operators, it is, however, still necessary to compute the relevant scheme matching factors perturbatively. We compute the latter at one loop, which, in particular, allows us to determine the complete set of next-to-leading (NLO) anomalous dimensions in our SF schemes.

An interesting byproduct of our computation is the possibility to study the systematic uncertainties related to the use of NLO perturbation theory in the computation of the RG running of four-quark operators to hadronic scales. This is a common feature of the phenomenological literature, and the question can be posed whether perturbative truncation effects can have an impact in physics analyses. The latter are studied in detail in our SF schemes, as well as in the ${\overline{\mathrm{MS}}}$ and RI-MOM schemes that have been studied in the literature. One of our main conclusions is that perturbative truncation effects in RG running can be argued to be significantly large. This makes a very strong case for a fully non-perturbative RG running program for these operators.

The structure of the paper is as follows. In Sect. 2 we provide a short review of the renormalization properties of the full basis of $\Delta F=2$ four-quark operators, stressing how considering it also allows one to obtain the anomalous dimensions of $\Delta F=1$ operators. We focus on the operators that appear in BSM physics, which exhibit scale-dependent mixing under renormalization, and discuss the definition of RGI operators in that case. In Sect. 3 we introduce our SF schemes, and explain the strategy to obtain NLO anomalous dimensions in the latter through a one-loop computation of the relevant four- and two-point correlation functions. Finally, in Sect. 4 we carry out a systematic study of the perturbative RG running in several schemes, and provide estimates of the resulting truncation uncertainty at scales in the few GeV range. In order to improve readability, several tables and figures are collected after the main text, and a many technical details are discussed in appendices.

2 Renormalization of four-quark operators

2.1 Mixing of four-quark operators under renormalization

The mixing under renormalization of four-quark operators that do not require subtraction of lower-dimensional operators has been determined in full generality in [1]. The absence of subtractions is elegantly implemented by using a formalism in which the operators are made of four different quark flavors; a complete set of Lorentz-invariant operators is

$$\begin{aligned} Q_1^\pm= & {} \mathcal{O}^\pm _\mathrm{VV+AA}, \quad \mathcal{Q}_1^\pm = \mathcal{O}^\pm _\mathrm{VA+AV},\nonumber \\ Q_2^\pm= & {} \mathcal{O}^\pm _\mathrm{VV-AA}, \quad \mathcal{Q}_2^\pm = \mathcal{O}^\pm _\mathrm{VA-AV},\nonumber \\ Q_3^\pm= & {} \mathcal{O}^\pm _\mathrm{SS-PP}, \quad \mathcal{Q}_3^\pm = \mathcal{O}^\pm _\mathrm{PS-SP},\\ Q_4^\pm= & {} \mathcal{O}^\pm _\mathrm{SS+PP}, \quad \mathcal{Q}_4^\pm = \mathcal{O}^\pm _\mathrm{SP+PS},\nonumber \\ Q_5^\pm= & {} \mathcal{O}^\pm _\mathrm{TT}, \quad \mathcal{Q}_5^\pm = \mathcal{O}^\pm _{\mathrm{T}\tilde{T}},\nonumber \end{aligned}$$

(2.1)

where

$$\begin{aligned} \mathcal{O}^\pm _{\Gamma _1\Gamma _2} = \frac{1}{2}\left[ (\bar{\psi }_1\Gamma _1\psi _2)(\bar{\psi }_3\Gamma _2\psi _4)\,\pm \, (\bar{\psi }_1\Gamma _1\psi _4)(\bar{\psi }_3\Gamma _2\psi _2) \right] ,\nonumber \\ \end{aligned}$$

(2.2)

$\mathcal{O}^\pm _{\Gamma _1\Gamma _2\pm \Gamma _2\Gamma _1}\equiv \mathcal{O}^\pm _{\Gamma _1\Gamma _2} \pm \mathcal{O}^\pm _{\Gamma _2\Gamma _1}$, and the labeling is adopted $\mathrm{V}\rightarrow \gamma _\mu$, $\mathrm{A}\rightarrow \gamma _\mu \gamma _5$, $\mathrm{S}\rightarrow \mathbf {1}$, $\mathrm{P}\rightarrow \gamma _5$, $\mathrm{T}\rightarrow \sigma _{\mu \nu }$, ${\tilde{\mathrm{T}}}\rightarrow {\scriptstyle {{1\over 2}}}\varepsilon _{\mu \nu \rho \tau }\sigma _{\rho \tau }$, with $\sigma _{\mu \nu }={\scriptstyle {{i\over 2}}}[\gamma _\mu ,\gamma _\nu ]$. In the above expression round parentheses indicate spin and color scalars, and subscripts are flavor labels. Note that operators $Q_k^\pm$ are parity-even, and $\mathcal{Q}_k^\pm$ are parity-odd.

It is important to stress that this framework is fairly general. For instance, with the assignments

$$\begin{aligned} \psi _1=\psi _3=s,\quad \psi _2=\psi _4=d \end{aligned}$$

(2.3)

the operators $Q_k^-$ vanish, while $Q_1^+$ enters the SM amplitude for $K^0$–$\bar{K}^0$ mixing, and $Q_{2,\ldots ,5}^+$ the contributions to the same amplitude from arbitrary extensions of the SM. Idem for $B_{(s)}^0$–$\bar{B}_{(s)}^0$ mixing with

$$\begin{aligned} \psi _1=\psi _3=b,\quad \psi _2=\psi _4=d/s. \end{aligned}$$

(2.4)

If one instead chooses the assignments

$$\begin{aligned}&\psi _1=s,\quad \psi _2=d,\quad \psi _3=\psi _4=u,c,\nonumber \\&\psi _1=s,\quad \psi _4=d,\quad \psi _2=\psi _3=u,c, \end{aligned}$$

(2.5)

the resulting $Q_1^\pm$ will be the operators in the SM $\Delta S=1$ effective weak Hamiltonian with an active charm quark, which, in the chiral limit, do not mix with lower-dimensional operators. By proceeding in this way, essentially all possible four-quark effective interactions with net flavor change can easily be seen to be comprised within our scope.

In the following we will assume a mass-independent renormalization scheme. Renormalized operators can be written as

(2.6)

(summations over l, m are implied), where the matrices $Z,\mathcal{Z}$ are scale dependent and reabsorb logarithmic divergences, while are (possible) finite subtraction coefficients that only depend on the bare coupling. They have the generic structure

$$\begin{aligned}&Z=\left( \begin{array}{ccccc} Z_{11} &{} 0 &{} 0 &{} 0 &{} 0 \\ 0 &{} Z_{22} &{} Z_{23} &{} 0 &{} 0 \\ 0 &{} Z_{32} &{} Z_{33} &{} 0 &{} 0 \\ 0 &{} 0 &{} 0 &{} Z_{44} &{} Z_{45} \\ 0 &{} 0 &{} 0 &{} Z_{54} &{} Z_{55} \end{array}\right) ,\nonumber \\&\Delta =\left( \begin{array}{ccccc} 0 &{} \Delta _{12} &{} \Delta _{13} &{} \Delta _{14} &{} \Delta _{15} \\ \Delta _{21} &{} 0 &{} 0 &{} \Delta _{24} &{} \Delta _{25} \\ \Delta _{31} &{} 0 &{} 0 &{} \Delta _{34} &{} \Delta _{35} \\ \Delta _{41} &{} \Delta _{42} &{} \Delta _{43} &{} 0 &{} 0 \\ \Delta _{51} &{} \Delta _{52} &{} \Delta _{53} &{} 0 &{} 0 \end{array}\right) , \end{aligned}$$

(2.7)

and similarly in the parity-odd sector. If chiral symmetry is preserved by the regularization, both $\Delta$ and vanish. The main result of [1] is that even when a lattice regularization that breaks chiral symmetry explicitly through the Wilson term is employed, due to the presence of residual discrete flavor symmetries. In particular, the left–left operators $\mathcal{Q}_\mathrm{VA+AV}^\pm$ that mediate Standard Model-allowed transitions renormalize multiplicatively, while operators that appear as effective interactions in extensions of the Standard Model do always mix.¹

Interestingly, in [1] some identities are derived that relate the renormalization matrices for $(\mathcal{Q}_2^+,\mathcal{Q}_3^+)$ and $(\mathcal{Q}_2^-,\mathcal{Q}_3^-)$ in RI-MOM schemes. In Appendix A we discuss the underlying symmetry structure in some more detail, and show how it can be used to derive constraints between matrices of anomalous dimensions in generic schemes.

2.2 Callan–Symanzik equations

Theory parameters and operators are renormalized at the renormalization scale $\mu$. The scale dependence of renormalized quantities is then governed by renormalization group evolution. We will consider QCD with $N_{\mathrm {\scriptstyle f}}$ quark flavors and $N$ colors. The Callan–Symanzik equations satisfied by the gauge coupling and quark masses are of the form

$$\begin{aligned} q\frac{\partial }{\partial q}\overline{g}(q)&= \beta (\overline{g}(q)),\end{aligned}$$

(2.8)

$$\begin{aligned} q\frac{\partial }{\partial q}\overline{m}_\mathrm{f}(q)&= \tau (\overline{g}(q))\overline{m}_\mathrm{f}(q), \end{aligned}$$

(2.9)

respectively, and satisfy the initial conditions

$$\begin{aligned} \overline{g}(\mu )&= g_\mathrm{R},\end{aligned}$$

(2.10)

$$\begin{aligned} \overline{m}_\mathrm{f}(\mu )&= m_{\mathrm{R} ,\mathrm f}, \end{aligned}$$

(2.11)

where $\mathrm{f}$ is a flavor label. Mass independence of the scheme is reflected in the fact that the beta function and mass anomalous dimension $\tau$ depend on the coupling and the number of flavors, but not on quark masses. Asymptotic perturbative expansions read

$$\begin{aligned} \beta (g)&\underset{g \sim 0}{\approx } -g^3(b_0+b_1g^2+\cdots ),\end{aligned}$$

(2.12)

$$\begin{aligned} \tau (g)&\underset{g \sim 0}{\approx } -g^2(d_0+d_1g^2+\cdots ). \end{aligned}$$

(2.13)

The universal coefficients of the perturbative beta function and mass anomalous dimension are

$$\begin{aligned} b_0= & {} \frac{1}{(4\pi )^2}\left[ \frac{11}{3}N-\frac{2}{3}N_{\mathrm {\scriptstyle f}}\right] ,\nonumber \\ b_1= & {} \frac{1}{(4\pi )^4}\left[ \frac{34}{3}N^2 -\left( \frac{13}{3}N-\frac{1}{N}\right) N_{\mathrm {\scriptstyle f}}\right] ,\\ d_0= & {} \frac{1}{(4\pi )^2}\,\frac{3(N^2-1)}{N}.\nonumber \end{aligned}$$

(2.14)

We will deal with Euclidean correlation functions of gauge-invariant composite operators. Without loss of generality, let us consider correlation functions of the form

$$\begin{aligned} G_k(x;y_1,\ldots ,y_n) = \langle O_k(x)\mathcal{O}_1(y_1)\cdots \mathcal{O}_n(y_n))\rangle , \end{aligned}$$

(2.15)

with $x \ne y_j~\forall j,~y_j \ne y_k~\forall j \ne k$, where $\{O_k\}$ is a set of operators that mix under renormalization, and where $\mathcal{O}_k$ are multiplicatively renormalizable operators.² Renormalized correlation functions satisfy a system of Callan–Symanzik equations obtained by imposing the independence of $G_k$ from the renormalization scale $\mu$, viz.

$$\begin{aligned} \mu \frac{\mathrm{d}}{\mathrm{d}\mu } \bar{G}_j = \sum _k\left[ \gamma _{jk}(g_\mathrm{R})-\sum _{l=1}^n\tilde{\gamma }_l(g_\mathrm{R})\right] \bar{G}_k, \end{aligned}$$

(2.16)

which, expanding the total derivative, leads to

$$\begin{aligned}&\Bigg \{ \mu \frac{\partial }{\partial \mu } + \beta (g_\mathrm{R})\frac{\partial }{\partial g_\mathrm{R}} + \beta _\lambda (g_\mathrm{R})\lambda \,\frac{\partial }{\partial \lambda } \nonumber \\&\quad + \sum _{\mathrm{f}=1}^{N_{\mathrm {\scriptstyle f}}}\tau (g_\mathrm{R})m_{\mathrm{R} ,\mathrm f}\frac{\partial }{\partial m_{\mathrm{R} ,\mathrm f}}\nonumber \\&\quad - \sum _{l=1}^n\tilde{\gamma }_l(g_\mathrm{R}) \Bigg \}\bar{G}_j = \sum _k\gamma _{jk}(g_\mathrm{R})\,\bar{G}_k , \end{aligned}$$

(2.17)

where $\gamma$ is a matrix of anomalous dimensions describing the mixing of $\{O_k\}$, and $\tilde{\gamma }_l$ is the anomalous dimension of $\mathcal{O}_l$. For completeness, we have included a term which takes into account the dependence on the gauge parameter $\lambda$ in covariant gauges; this term is absent in schemes like ${\overline{\mathrm{MS}}}$ (irrespective of the regularization prescription) or the SF schemes we will introduce, but it is present in the RI schemes we will also be concerned with later. The RG function $\beta _\lambda$ is given by

$$\begin{aligned} q\frac{\partial }{\partial q}\lambda (q) = \beta _\lambda (\overline{g}(q))\lambda (q), \end{aligned}$$

(2.18)

and its perturbative expansion has the form

$$\begin{aligned} \beta _\lambda (g) = -g^2(b_0^\lambda + b_1^\lambda g^2 + \cdots ), \end{aligned}$$

(2.19)

where the universal coefficient is given by

$$\begin{aligned} b_0^\lambda = \frac{1}{(4\pi )^2}\left[ \left( \lambda -\frac{13}{3}\right) N+\frac{4}{3}N_{\mathrm {\scriptstyle f}}\right] . \end{aligned}$$

(2.20)

In the Landau gauge ($\lambda =0$) the term with $\beta _\lambda$ always vanishes. From now on, in order to avoid unnecessary complications, we will assume that whenever RI anomalous dimensions are employed they will be in Landau gauge, and consequently drop terms with $\beta _\lambda$ in all equations.

From now on, in order to simplify the notation we will use the shorthand notation

$$\begin{aligned} q\frac{\partial }{\partial q}\overline{O}_j(q) = \sum _k\gamma _{jk}(\overline{g}(q))\overline{O}_k(q) \end{aligned}$$

(2.21)

for the Callan–Symanzik equation satisfied by the insertion of a composite operator in a renormalized, on-shell correlation function (i.e. Eq. (2.21) is to be interpreted in the sense provided by Eq. (2.17)). The corresponding initial condition can be written as

$$\begin{aligned} \overline{O}_k(\mu ) = {O_\mathrm{R}}_{,k}, \end{aligned}$$

(2.22)

and the perturbative expansion of the anomalous dimension matrix $\gamma$ as

$$\begin{aligned} \gamma (g) \underset{g \sim 0}{\approx } -g^2(\gamma _0+\gamma _1g^2+\cdots ). \end{aligned}$$

(2.23)

The universal, one-loop coefficients of the anomalous dimension matrix for four-fermion operators were first computed in [5,6,7]. With our notational conventions the non-zero entries read

$$\begin{aligned} \gamma ^{+,(0)}_{11}= & {} \left( 6-\frac{6}{N} \right) (4\pi )^{-2}, \nonumber \\ \gamma ^{-,(0)}_{11}= & {} \left( -6-\frac{6}{N}\right) (4\pi )^{-2} ,\nonumber \\ \gamma ^{+,(0)}_{22}= & {} \left( \frac{6}{N}\right) (4\pi )^{-2} , \quad \gamma ^{-,(0)}_{22} =\left( \frac{6}{N}\right) (4\pi )^{-2} ,\nonumber \\ \gamma ^{+,(0)}_{23}= & {} 12 (4\pi )^{-2},\quad \gamma ^{-,(0)}_{23} = -12(4\pi )^{-2}, \nonumber \\ \gamma ^{+,(0)}_{33}= & {} \left( -6N+\frac{6}{N}\right) (4\pi )^{-2} \nonumber \\ \gamma ^{-,(0)}_{33}= & {} \left( -6N+\frac{6}{N}\right) (4\pi )^{-2} ,\nonumber \\ \gamma ^{+,(0)}_{44}= & {} \left( 6-6N+\frac{6}{N}\right) (4\pi )^{-2} ,\nonumber \\ \gamma ^{-,(0)}_{44}= & {} \left( -6-6N+\frac{6}{N}\right) (4\pi )^{-2} ,\nonumber \\ \gamma ^{+,(0)}_{45}= & {} \left( \frac{1}{2}-\frac{1}{N}\right) (4\pi )^{-2},\\ \gamma ^{-,(0)}_{45}= & {} \left( -\frac{1}{2}-\frac{1}{N}\right) (4\pi )^{-2},\nonumber \\ \gamma ^{+,(0)}_{54}= & {} \left( -24-\frac{48}{N}\right) (4\pi )^{-2} , \nonumber \\ \gamma ^{-,(0)}_{54}= & {} \left( 24-\frac{48}{N}\right) (4\pi )^{-2} ,\nonumber \\ \gamma ^{+,(0)}_{55}= & {} \left( 6+2N-\frac{2}{N}\right) (4\pi )^{-2} , \nonumber \\ \gamma ^{-,(0)}_{55}= & {} \left( -6+2N-\frac{2}{N}\right) (4\pi )^{-2} .\nonumber \end{aligned}$$

(2.24)

2.3 Formal solution of the RG equation

Let us now consider the solution to Eq. (2.21). For that purpose we start by introducing the (matricial) renormalization group evolution operator $U(\mu _2,\mu _1)$ that evolves renormalized operators between the scales³ $\mu _1$ and $\mu _2<\mu _1$,

$$\begin{aligned} \overline{O}_i(\mu _2) = U_{ij}(\mu _2,\mu _1) \overline{O}_j(\mu _1). \end{aligned}$$

(2.25)

By substituting into Eq. (2.21) one has the equation for $U(\mu _2,\mu _1)$

$$\begin{aligned} \mu _2\,\frac{\partial }{\partial \mu _2}\,U(\mu _2,\mu _1) = \gamma [\overline{g}(\mu _2)]U(\mu _2,\mu _1) \end{aligned}$$

(2.26)

(n.b. the matrix product on the r.h.s.) with initial condition $U(\mu _1,\mu _1)=\mathbf {1}$. Following a standard procedure, this differential equation for U can be converted into a Volterra-type integral equation and solved iteratively, viz.

$$\begin{aligned} U(\mu _2,\mu _1) = \mathrm{T}\exp \left\{ \int _{\overline{g}(\mu _1)}^{\overline{g}(\mu _2)}\mathrm{d}g\,\frac{1}{\beta (g)}\,\gamma (g) \right\} , \end{aligned}$$

(2.27)

where as usual the notation $\mathrm{T}\exp$ refers to a definition in terms of the Taylor expansion of the exponential function with “powers” of the integral involving argument-ordered integrands – explicitly, for a generic matrix function M, one has

$$\begin{aligned}&\mathrm{T}\exp \left\{ \int _{x_-}^{x_+}\mathrm{d}x\,M(x)\right\} \equiv \mathbf {1} + \int _{x_-}^{x_+}\mathrm{d}x\,M(x)\nonumber \\&\qquad + \int _{x_-}^{x_+}\mathrm{d}x_1\,M(x_1)\int _{x_-}^{x_1}\mathrm{d}x_2\,M(x_2)\nonumber \\&\qquad + \int _{x_-}^{x_+}\mathrm{d}x_1\,M(x_1)\int _{x_-}^{x_1}\mathrm{d}x_2\,M(x_2)\int _{x_-}^{x_2}\mathrm{d}x_3\,M(x_3)\nonumber \\&\qquad + \cdots = \mathbf {1} + \int _{x_-}^{x_+}\mathrm{d}x\,M(x)\nonumber \\&\qquad + \frac{1}{2!}\int _{x_-}^{x_+}\mathrm{d}x_1\int _{x_-}^{x_+}\mathrm{d}x_2\,\Big \{\theta (x_1-x_2)M(x_1)M(x_2)\nonumber \\&\qquad +\theta (x_2-x_1)M(x_2)M(x_1)\Big \}+ \cdots \end{aligned}$$

(2.28)

2.4 RGI in the absence of mixing

Let us now consider an operator O that renormalizes multiplicatively. In that case, both $\gamma$ and U are scalar objects, and Eq. (2.25) can be manipulated as

$$\begin{aligned}&\overline{O}(\mu _2) = \exp \left\{ \int _{\overline{g}(\mu _1)}^{\overline{g}(\mu _2)}\mathrm{d}g\,\frac{\gamma (g)}{\beta (g)} \right\} \overline{O}(\mu _1)\nonumber \\&\quad =\exp \left\{ \int _{\overline{g}(\mu _1)}^{\overline{g}(\mu _2)}\mathrm{d}g\,\frac{\gamma _0}{b_0g} \right\} \nonumber \\&\qquad \times \, \exp \left\{ \int _{\overline{g}(\mu _1)}^{\overline{g}(\mu _2)}\mathrm{d}g\,\left[ \frac{\gamma (g)}{\beta (g)}-\frac{\gamma _0}{b_0g}\right] \right\} \overline{O}(\mu _1)\nonumber \\&\quad =\left[ \frac{\overline{g}^2(\mu _2)}{\overline{g}^2(\mu _1)}\right] ^{\frac{\gamma _0}{2b_0}} \exp \left\{ \int _{\overline{g}(\mu _1)}^{\overline{g}(\mu _2)}\mathrm{d}g\,\left[ \frac{\gamma (g)}{\beta (g)}-\frac{\gamma _0}{b_0g}\right] \right\} \overline{O}(\mu _1), \end{aligned}$$

(2.29)

yielding the identity⁴

$$\begin{aligned}&\left[ \frac{\overline{g}^2(\mu _2)}{4\pi }\right] ^{-\frac{\gamma _0}{2b_0}}\overline{O}(\mu _2) =\left[ \frac{\overline{g}^2(\mu _1)}{4\pi }\right] ^{-\frac{\gamma _0}{2b_0}}\nonumber \\&\quad \times \exp \left\{ \int _{\overline{g}(\mu _1)}^{\overline{g}(\mu _2)}\mathrm{d}g\,\left[ \frac{\gamma (g)}{\beta (g)}-\frac{\gamma _0}{b_0g}\right] \right\} \overline{O}(\mu _1). \end{aligned}$$

(2.30)

The advantage of having rewritten Eq. (2.25) in this way is that now the integral in the exponential is finite as either integration limit is taken to zero; in particular, the r.h.s. is well defined when $\mu _2\rightarrow \infty ~\Leftrightarrow ~\overline{g}(\mu _2)\rightarrow 0$, and therefore so is the l.h.s. Thus, we define the RGI operator insertion as

$$\begin{aligned} \hat{O} \equiv \lim _{\mu \rightarrow \infty }\left[ \frac{\overline{g}^2(\mu )}{4\pi }\right] ^{-\frac{\gamma _0}{2b_0}} \overline{O}(\mu ), \end{aligned}$$

(2.31)

upon which we have an explicit expression to retrieve the RGI operator from the renormalized one at any value of the renormalization scale $\mu$, provided the anomalous dimension and the beta function are known for scales $\ge \mu$,

$$\begin{aligned} \hat{O} {=} \left[ \frac{\overline{g}^2(\mu )}{4\pi }\right] ^{-\frac{\gamma _0}{2b_0}} \exp \left\{ {-}\int _{0}^{\overline{g}(\mu )}\mathrm{d}g\,\left[ \frac{\gamma (g)}{\beta (g)}{-}\frac{\gamma _0}{b_0g}\right] \right\} \overline{O}(\mu ).\nonumber \\ \end{aligned}$$

(2.32)

Starting from the latter equation, it is easy to check explicitly that $\hat{O}$ is invariant under a change of renormalization scheme.

Note that the crucial step in the manipulation has been to add and subtract the term $\frac{\gamma _0}{b_0g}$ in the integral that defines the RG evolution operator, which allows one to obtain a quantity that is UV-finite by removing the logarithmic divergence induced at small g by the perturbative behavior $\gamma (g)/\beta (g) \sim 1/g$. When $\gamma$ is a matrix of anomalous dimensions this step becomes non-trivial, since in general $[\gamma (g),\gamma _0]\ne 0$; the derivation has thus to be performed somewhat more carefully.

2.5 RGI in the presence of mixing

Let us start by studying the UV behavior of the matricial RG evolution operator in Eq. (2.25), using its actual definition in Eq. (2.28). For that purpose, we first observe that by taking the leading-order approximation for $\gamma (g)/\beta (g)$ the T-exponential becomes a standard exponential, since $[\gamma _0 g_1^2,\gamma _0 g_2^2]=0~\forall g_1,g_2$. One can then perform the integral trivially and write

$$\begin{aligned} U(\mu _2,\mu _1) \underset{\mathrm{LO}}{=} \left[ \frac{\overline{g}^2(\mu _2)}{\overline{g}^2(\mu _1)}\right] ^{\frac{\gamma _0}{2b_0}} \equiv U_\mathrm{\scriptscriptstyle LO}(\mu _2,\mu _1). \end{aligned}$$

(2.33)

When next-to-leading order corrections are included the T-exponential becomes non-trivial. In order to make contact with the literature (see e.g. [5, 8]), we write⁵

$$\begin{aligned} U(\mu _2,\mu _1) \equiv \left[ W(\mu _2)\right] ^{-1}\, U_\mathrm{\scriptscriptstyle LO}(\mu _2,\mu _1) W(\mu _1). \end{aligned}$$

(2.34)

Upon inserting Eqs. (2.35) in (2.26) we obtain for W the RG equation

$$\begin{aligned} \mu \frac{\partial }{\partial \mu }W(\mu )= & {} -W(\mu )\gamma (\overline{g}(\mu ))+\beta (\overline{g}(\mu )) \frac{\gamma _0}{b_0\overline{g}(\mu )}W(\mu ) \nonumber \\= & {} [\gamma (\overline{g}(\mu )),W(\mu )] \nonumber \\&-\, \beta (\overline{g}(\mu ))\left( \frac{\gamma (\overline{g}(\mu ))}{\beta (\overline{g}(\mu ))}-\frac{\gamma _0}{b_0\overline{g}(\mu )} \right) W(\mu ).\nonumber \\ \end{aligned}$$

(2.36)

The matrix W can be interpreted as the piece of the evolution operator containing contributions beyond the leading perturbative order. It is easy to check by expanding perturbatively (see below) that W is regular in the UV, and that all the logarithmic divergences in the evolution operator are contained in $U_\mathrm{\scriptscriptstyle LO}$; in particular,

$$\begin{aligned} W(\mu )\underset{\mu \rightarrow \infty }{=}\mathbf {1}. \end{aligned}$$

(2.37)

Note also that in the absence of mixing Eq. (2.36) can be solved explicitly to get (using $W(0)=1$)

$$\begin{aligned} W(\mu )\underset{\mathrm{no~mixing}}{=}\exp \left\{ -\int _{0}^{\overline{g}(\mu )}\mathrm{d}g\,\left[ \frac{\gamma (g)}{\beta (g)}-\frac{\gamma _0}{b_0g} \right] \right\} . \end{aligned}$$

(2.38)

Now it is easy, by analogy with the non-mixing case, to define RGI operators. We rewrite Eq. (2.25) as

$$\begin{aligned}&\left[ \frac{\overline{g}^2(\mu _2)}{4\pi }\right] ^{-\frac{\gamma _0}{2b_0}} W(\mu _2)\overline{O}(\mu _2) \nonumber \\&\quad =\left[ \frac{\overline{g}^2(\mu _1)}{4\pi }\right] ^{-\frac{\gamma _0}{2b_0}} W(\mu _1)\overline{O}(\mu _1), \end{aligned}$$

(2.39)

where $\overline{O}$ is a vector of renormalized operators on which the RG evolution matrix acts, cf. Eq. (2.25). The l.h.s. (resp. r.h.s.) is obviously finite for $\mu _1\rightarrow \infty$ (resp. $\mu _2\rightarrow \infty$), which implies that the vector of RGI operators can be obtained:

$$\begin{aligned} \hat{O} = \left[ \frac{\overline{g}^2(\mu )}{4\pi }\right] ^{-\frac{\gamma _0}{2b_0}}W(\mu )\overline{O}(\mu ) \equiv \tilde{U}(\mu )\overline{O}(\mu ). \end{aligned}$$

(2.40)

When there is no mixing, the use of Eq. (2.38) immediately brings back Eq. (2.32).

2.6 Perturbative expansion of RG evolution functions

By expanding Eq. (2.36) perturbatively, with⁶

$$\begin{aligned}&W(\mu ) \approx \mathbf {1} + \overline{g}^2(\mu )J_1 + \overline{g}^4(\mu )J_2 \nonumber \\&\qquad \qquad +\, \overline{g}^6(\mu )J_3 + \overline{g}^8(\mu )J_4 + \cdots \end{aligned}$$

(2.41)

we find for the first four orders in the expansion the conditions

$$\begin{aligned} 2b_0J_1-\left[ \gamma _0,J_1\right]&= \frac{b_1}{b_0}\gamma _0-\gamma _1,\end{aligned}$$

(2.42)

$$\begin{aligned} 4b_0J_2-\left[ \gamma _0,J_2\right]&= J_1\left( \frac{b_1}{b_0}\gamma _0-\gamma _1\right) \nonumber \\&\quad +\left( \frac{b_2}{b_0}-\frac{b_1^2}{b_0^2}\right) \gamma _0 +\frac{b_1}{b_0}\gamma _1-\gamma _2,\end{aligned}$$

(2.43)

$$\begin{aligned} 6b_0J_3-\left[ \gamma _0,J_3\right]&= J_2\left( \frac{b_1}{b_0}\gamma _0-\gamma _1\right) \nonumber \\&\quad +J_1\left\{ \left( \frac{b_2}{b_0}-\frac{b_1^2}{b_0^2}\right) \gamma _0 +\frac{b_1}{b_0}\gamma _1-\gamma _2\right\} \nonumber \\&\quad +\left( \frac{b_3}{b_0}-\frac{2b_2b_1}{b_0^2}+\frac{b_1^3}{b_0^3}\right) \gamma _0\nonumber \\&\quad +\left( \frac{b_2}{b_0}-\frac{b_1^2}{b_0^2}\right) \gamma _1 +\frac{b_1}{b_0}\gamma _2-\gamma _3,\end{aligned}$$

(2.44)

$$\begin{aligned} 8b_0J_4-\left[ \gamma _0,J_4\right]&= J_3\left( \frac{b_1}{b_0}\gamma _0-\gamma _1\right) \nonumber \\&\quad +J_2\left\{ \left( \frac{b_2}{b_0}-\frac{b_1^2}{b_0^2}\right) \gamma _0 +\frac{b_1}{b_0}\gamma _1-\gamma _2\right\} \nonumber \\&\quad +J_1\left\{ \left( \frac{b_3}{b_0}-\frac{2b_2b_1}{b_0^2} +\frac{b_1^3}{b_0^3}\right) \gamma _0\right. \nonumber \\&\quad \left. +\left( \frac{b_2}{b_0} -\frac{b_1^2}{b_0^2}\right) \gamma _1+\frac{b_1}{b_0}\gamma _2 -\gamma _3\right\} \nonumber \\&\quad +\left( -\frac{b_1^4}{b_0^4} + 3\frac{b_1^2b_2}{b_0^3} -\frac{b_2^2}{b_0^2} - 2\frac{b_1b_3}{b_0^2} + \frac{b_4}{b_0}\right) \gamma _0\nonumber \\&\quad +\left( \frac{b_3}{b_0}-\frac{2b_2b_1}{b_0^2} +\frac{b_1^3}{b_0^3}\right) \gamma _1\nonumber \\&\quad +\left( \frac{b_2}{b_0}-\frac{b_1^2}{b_0^2}\right) \gamma _2 +\frac{b_1}{b_0}\gamma _3-\gamma _4. \end{aligned}$$

(2.45)

Modulo sign and normalization conventions (involving powers of $4\pi$ related to expanding in $\overline{g}^2$ rather than $\alpha /(4\pi )$), and the dependence on gauge fixing (which does not apply to our context), Eq. (2.42) coincides with Eq. (24) of [5]. All four equations, as well as those for higher orders, can easily be solved to obtain $J_n$ for given values the coefficients in the perturbative expansion of $\gamma$. The LO, NLO, and NNLO and NNNLO matching for the RGI operators is thus obtained from Eq. (2.40) by using the expansion in powers of $\overline{g}^2$ in Eq. (2.41) up to zeroth, first, second, and third order, respectively.

3 Anomalous dimensions in SF schemes

3.1 Changes of renormalization scheme

Let us now consider a change to a different mass-independent renormalization scheme, indicated by primes. The relation between renormalized quantities in either scheme amounts to finite renormalizations of the form

$$\begin{aligned} g_\mathrm{R}'&= \sqrt{\mathcal{X}_\mathrm{g}(g_\mathrm{R})}\,g_\mathrm{R},\end{aligned}$$

(3.1)

$$\begin{aligned} m_{\mathrm{R} ,\mathrm f}'&= \mathcal{X}_\mathrm{m}(g_\mathrm{R})\,m_{\mathrm{R} ,\mathrm f},\end{aligned}$$

(3.2)

$$\begin{aligned} {O_\mathrm{R}}_{,j}'&= (\mathcal{X}_{O})_{jk}\,{O_\mathrm{R}}_{,k}. \end{aligned}$$

(3.3)

The scheme-change factors $\mathcal{X}$ can be expanded perturbatively as

$$\begin{aligned} \mathcal{X}(g) \underset{g \sim 0}{\approx } 1 + \mathcal{X}^{(1)} g^2 + \cdots . \end{aligned}$$

(3.4)

By substituting Eqs. (3.1–3.3) into the corresponding Callan–Symanzik equations, the relation between the RG evolution functions in different schemes is found,

$$\begin{aligned} \beta '(g')&= \left[ \beta (g)\frac{\partial g'}{\partial g}\right] _{g=g(g')},\end{aligned}$$

(3.5)

$$\begin{aligned} \tau '(g')&= \left[ \tau (g)+\beta (g)\frac{\partial }{\partial g}\ln \mathcal{X}_\mathrm{m}(g)\right] _{g=g(g')},\end{aligned}$$

(3.6)

$$\begin{aligned} \gamma '(g')&= \left[ \gamma (g)+\beta (g)\frac{\partial }{\partial g}\ln \mathcal{X}_O(g)\right] _{g=g(g')}. \end{aligned}$$

(3.7)

One can then plug in perturbative expansions and obtain explicit formulae relating coefficients in different schemes. In particular, it is found that $b_0,b_1$ are scheme-independent, and the same applies to $d_0$ and $\gamma ^{(0)}$. The relation between next-to-leading order coefficients for quark masses and operator anomalous dimensions are given by

$$\begin{aligned} d_1'&= d_1 + 2b_0\mathcal{X}_\mathrm{m}^{(1)}-d_0\mathcal{X}_\mathrm{g}^{(1)},\end{aligned}$$

(3.8)

$$\begin{aligned} {\gamma '}^{(1)}&= \gamma ^{(1)} + \left[ \mathcal{X}_O^{(1)},\gamma ^{(0)}\right] \nonumber \\&\quad + 2b_0\mathcal{X}_O^{(1)} + b_0^\lambda \frac{\partial }{\partial \lambda }\mathcal{X}_O^{(1)} - \gamma ^{(0)}\mathcal{X}_\mathrm{g}^{(1)}. \end{aligned}$$

(3.9)

Therefore, if the anomalous dimension is known at two loops in some scheme, in order to obtain the same result in a different scheme it is sufficient to compute the one-loop relation between them.

3.2 Strategy for the computation of NLO anomalous dimensions in SF schemes

Equation (3.9) will be the key ingredient for our computation of anomalous dimensions to two loops in SF schemes, using as starting point known two-loop results in ${\overline{\mathrm{MS}}}$ or RI schemes. Indeed, our strategy will be to compute the one-loop matching coefficient between the SF schemes that we will introduce presently, and the continuum schemes where $\gamma ^{(1)}$ is known. $\gamma ^{(1);{\overline{\mathrm{MS}}}}$ can be found in [8,9,10], while $\gamma ^{(1);\mathrm RI}$ can be computed from both [5, 8]; we gather them in Appendix B.

One practical problem arises due to the dependence of the scheme definition in the continuum on the regulator employed (usually some form of dimensional regularization). This implies that one-loop computation in SF schemes needed to obtain the matching coefficient should be carried out using the same regulator as in the continuum scheme. However, the lattice is the only regularization currently available for the SF. As a consequence, it is necessary to employ a third, intermediate reference scheme, which we will dub “lat”, where the ${\overline{\mathrm{MS}}}$ or RI prescription is applied to the lattice-regulated theory. One can then proceed in two steps:

• (i) Compute the matching coefficient $[\mathcal{X}_O^{(1)}]_\mathrm{SF;lat}$ between SF and lat schemes. As we will see later, the latter is retrieved by computing SF renormalization constants at one loop.

• (ii) Retrieve the one-loop matching coefficients between the lattice- and dimensionally-regulated versions of the continuum scheme “cont” (i.e. ${\overline{\mathrm{MS}}}$ or RI), $[\mathcal{X}_O^{(1)}]_\mathrm{cont;lat}$, and obtain the matching coefficient that enters Eq. (3.9) as

  $$\begin{aligned} \left[ \mathcal{X}_O^{(1)}\right] _\mathrm{SF;cont} = \left[ \mathcal{X}_O^{(1)}\right] _\mathrm{SF;lat}-\left[ \mathcal{X}_O^{(1)}\right] _\mathrm{cont;lat}. \end{aligned}$$

  (3.10)

The one-loop matching coefficients $[\mathcal{X}_O^{(1)}]_\mathrm{cont;lat}$ that we will need can be extracted from the literature. For the RI-MOM scheme they can be found in [11], while for the ${\overline{\mathrm{MS}}}$ scheme they can be extracted from [12,13,14]).⁷ We gather the RI-MOM results in Landau gauge in Appendix D. $\chi ^{(1)}_\mathrm{g}$ can be found in [15].

3.3 SF renormalization conditions

We now consider the problem of specifying suitable renormalization conditions on four-quark operators, using the Schrödinger Functional formalism. The latter [16], initially developed to produce a precise determination of the running coupling [17,18,19,20,21,22], has been extended along the years to various other phenomenological contexts, like e.g. quark masses [23,24,25] or heavy-light currents relevant for B-physics, among others [26, 27]. In the context of four-quark operators, the first applications involved the multiplicatively renormalizable operators $\mathcal{Q}_1^\pm$ of Eq. (2.1) (which, as explained above, enter Standard Model effective Hamiltonians for $\Delta F=1$ and $\Delta F=2$ processes) [28,29,30,31], as well as generic $\Delta B=2$ operators in the static limit [31, 32]. The latter studies are extended in this paper to cover the full set of relativistic operators. It is important to stress that, while these schemes will be ultimately employed in the context of a non-perturbative lattice computation of renormalization constants and anomalous dimensions, the definition of the schemes is fully independent of any choice of regulator.

We use the standard SF setup as described in [33], where the reader is referred for full details including unexplained notation. We will work on boxes with spatial extent L and time extent T; in practice, $T=L$ will always be set. Source fields are made up of boundary quarks and antiquarks,

$$\begin{aligned} \mathcal{O}_{\alpha \beta }[\Gamma ]&= a^6\sum _{\mathbf {y},\mathbf {z}}\bar{\zeta }_\alpha (\mathbf {y})\Gamma \zeta _\beta (\mathbf {z}),\end{aligned}$$

(3.11)

$$\begin{aligned} \mathcal{O}'_{\alpha \beta }[\Gamma ]&= a^6\sum _{\mathbf {y},\mathbf {z}}\bar{\zeta }'_\alpha (\mathbf {y})\Gamma \zeta '_\beta (\mathbf {z}), \end{aligned}$$

(3.12)

where $\alpha ,\beta$ are flavor indices, unprimed (primed) fields live at the $x_0=0$ ($x_0=T$) boundary, and $\Gamma$ is a spin matrix that must anticommute with $\gamma _0$, so that the boundary fermion field does not vanish. This is a consequence of the structure of the conditions imposed on boundary fields,

$$\begin{aligned} \zeta (\mathbf {x})={\scriptstyle {{1\over 2}}}(\mathbf {1}-\gamma _0)\zeta (\mathbf {x}),\quad \bar{\zeta }(\mathbf {x})=\bar{\zeta }(\mathbf {x}){\scriptstyle {{1\over 2}}}(\mathbf {1}+\gamma _0), \end{aligned}$$

(3.13)

and similarly for primed fields. The resulting limitations on the possible Dirac structures for these operators imply e.g. that it is not possible to use scalar bilinear operators, unless non-vanishing angular momenta are introduced. This can, however, be expected to lead to poor numerical behavior; thus, following our earlier studies [28, 29, 31, 32], we will work with zero-momentum bilinears and combine them suitably to produce the desired quantum numbers.

Renormalization conditions will be imposed in the massless theory, in order to obtain a mass-independent scheme by construction. They will furthermore be imposed on parity-odd four-quark operators, since working in the parity-even sector would entail dealing with the extra mixing due to explicit chiral symmetry breaking with Wilson fermions, cf. Eq. (2.7). In order to obtain non-vanishing SF correlation functions, we will then need a product of source operators with overall negative parity; taking into account the above observation about boundary fields, and the need to saturate flavor indices, the minimal structure involves three boundary bilinear operators and the introduction of an extra, “spectator” flavor (labeled number 5, keeping with the notation in Eq. (2.2)). We thus end up with correlation functions of the generic form

$$\begin{aligned} F_{k;s}^\pm (x_0)&= \langle \mathcal{Q}_k^\pm (x) \mathcal{S}_s \rangle ,\end{aligned}$$

(3.14)

$$\begin{aligned} G_{k;s}^\pm (x_0)&= \eta _k\langle \mathcal{Q}_k^\pm (x) \mathcal{S}'_s \rangle , \end{aligned}$$

(3.15)

where $\mathcal{S}_s$ is one of the five source operators

$$\begin{aligned} \mathcal{S}_1&= \mathcal{W}[\gamma _5,\gamma _5,\gamma _5],\end{aligned}$$

(3.16)

$$\begin{aligned} \mathcal{S}_2&= \frac{1}{6}\sum _{k,l,m=1}^3\epsilon _{klm} \mathcal{W}[\gamma _k,\gamma _l,\gamma _m],\end{aligned}$$

(3.17)

$$\begin{aligned} \mathcal{S}_3&= \frac{1}{3}\sum _{k=1}^3 \mathcal{W}[\gamma _5,\gamma _k,\gamma _k],\end{aligned}$$

(3.18)

$$\begin{aligned} \mathcal{S}_4&= \frac{1}{3}\sum _{k=1}^3 \mathcal{W}[\gamma _k,\gamma _5,\gamma _k],\end{aligned}$$

(3.19)

$$\begin{aligned} \mathcal{S}_5&= \frac{1}{3}\sum _{k=1}^3 \mathcal{W}[\gamma _k,\gamma _k,\gamma _5] \end{aligned}$$

(3.20)

with

$$\begin{aligned} \mathcal{W}[\Gamma _1,\Gamma _2,\Gamma _3] = L^{-3}\mathcal{O}'_{21}[\Gamma _1]\mathcal{O}'_{45}[\Gamma _2]\mathcal{O}_{53}[\Gamma _3], \end{aligned}$$

(3.21)

and similarly for $\mathcal{S}'_s$. The constant $\eta _k$ is a sign that ensures $F_{k;s}^\pm (x_0)=G_{k;s}^\pm (x_0)$ for all possible indices; it is easy to check that $\eta _2=-1,~\eta _{s \ne 2}=+1$.⁸ We will also use the two-point functions of boundary sources,

$$\begin{aligned} f_1&= -\,\frac{1}{2L^6}\langle \mathcal{O}'_{21}[\gamma _5]\mathcal{O}_{12} [\gamma _5]\rangle ,\end{aligned}$$

(3.22)

$$\begin{aligned} k_1&= -\,\frac{1}{6L^6}\sum _{k=1}^3\langle \mathcal{O}'_{21} [\gamma _k]\mathcal{O}_{12}[\gamma _k]\rangle . \end{aligned}$$

(3.23)

Finally, we define the ratios

$$\begin{aligned} \mathcal{A}_{k;s,\alpha }^\pm = \frac{F_{1;s}^\pm (T/2)}{f_1^{\scriptscriptstyle {\frac{3}{2}} -\alpha }k_1^\alpha }, \end{aligned}$$

(3.24)

where $\alpha$ is an arbitrary real parameter. The structure of $F_{k;s}^\pm$ and $f_1,k_1$ is illustrated in Fig. 1.

Fig. 1

Feynman diagrams for the four-quark correlation functions $F_{k;s}^\pm$ and the boundary-to-boundary correlators $f_1,k_1$ at tree level. Euclidean time goes from left to right. The double blob indicates the four-quark operator insertion, and dashed lines indicate the explicit time-like link variable involved in boundary-to-boundary quark propagators

We then proceed to impose renormalization conditions at bare coupling $g_0$ and scale $\mu =1/L$ by generalizing the condition introduced in [28, 29] for the renormalizable multiplicative operators $\mathcal{Q}_1^\pm$: the latter reads

$$\begin{aligned} Z_{1;s,\alpha }^\pm \,\mathcal{A}_{k;s,\alpha }^\pm = \left. \mathcal{A}_{k;s,\alpha }^\pm \right| _{g_0^2=0}, \end{aligned}$$

(3.25)

while, for operators that mix in doublets, we impose⁹

$$\begin{aligned}&\left( \begin{array}{cc} Z_{22;s_1,s_2,\alpha }^\pm &{} Z_{23;s_1,s_2,\alpha }^\pm \\ Z_{32;s_1,s_2,\alpha }^\pm &{} Z_{33;s_1,s_2,\alpha }^\pm \\ \end{array}\right) \left( \begin{array}{cc} \mathcal{A}_{2;s_1,\alpha }^\pm &{} \mathcal{A}_{2;s_2,\alpha }^\pm \\ \mathcal{A}_{3;s_1,\alpha }^\pm &{} \mathcal{A}_{3;s_2,\alpha }^\pm \\ \end{array}\right) \nonumber \\&\qquad = \left( \begin{array}{cc} \mathcal{A}_{2;s_1,\alpha }^\pm &{} \mathcal{A}_{2;s_2,\alpha }^\pm \\ \mathcal{A}_{3;s_1,\alpha }^\pm &{} \mathcal{A}_{3;s_2,\alpha }^\pm \\ \end{array}\right) _{g_0^2=0} , \end{aligned}$$

(3.26)

and similarly for $\mathcal{Q}_{4,5}^\pm$. The products of boundary-to-boundary correlators in the denominator of Eq. (3.24) cancel the renormalization of the boundary operators in $F_{k;s}^\pm$, and therefore $Z_{k;s,\alpha }^\pm$ only contains anomalous dimensions of four-fermion operators. Following [23, 28, 31], conditions are imposed on renormalization functions evaluated at $x_0=T/2$, and the phase that parameterizes spatial boundary conditions on fermion fields is fixed to $\theta =0.5$. Together with the $L=T$ geometry of our finite box, this fixes the renormalization scheme completely, up to the choice of boundary source, indicated by the index s, and the parameter $\alpha$. The latter can in principle take any value, but we will restrict our study to the choices $\alpha =0,1,3/2$.

One still has to check that renormalization conditions are well defined at tree level. While this is straightforward for Eq. (3.25), it is not so for Eq. (3.26): it is still possible that the matrix of ratios $\mathcal{A}$ has zero determinant at tree level, rendering the system of equations for the matrix of renormalization constants ill-conditioned. This is indeed the obvious case for $s_1 = s_2$, but the determinant turns out to be zero also for other non-trivial choices $s_1 \ne s_2$. In practice, out of the ten possible schemes one is only left with six, viz.¹⁰

$$\begin{aligned} (s_1,s_2) \in \{(1,2),(1,4),(1,5),(2,3),(3,4),(3,5)\}. \end{aligned}$$

(3.27)

It has to be stressed that this property is independent of the choice of $\theta$ and $\alpha$. Thus, we are left with a total of 15 schemes for $\mathcal{Q}_1^\pm$, and 18 for each of the pairs $(\mathcal{Q}_2^\pm ,\mathcal{Q}_3^\pm )$ and $(\mathcal{Q}_4^\pm ,\mathcal{Q}_5^\pm )$.

3.4 One-loop results in the SF

Let us now carry out a perturbative computation of the SF renormalization matrices introduced above, using a lattice regulator. For any of the correlation functions discussed in Sect. 3, the perturbative expansion reads

$$\begin{aligned} X=X^{(0)} + g_0^2\left[ X^{(1)}+m_\mathrm{cr}^{(1)}\frac{\partial X^{(0)}}{\partial m_0} \right] +\mathcal{O}(g_0^4), \end{aligned}$$

(3.28)

where X is one of $F_{k;s}^\pm (x_0),~f_1,~k_1$, or some combination thereof; where $m_0$ is the bare quark mass; and $m_\mathrm{cr}^{(1)}$ the one-loop coefficient in the perturbative expansion of the critical mass. The derivative term in the square bracket is needed to set the correlation function X to zero renormalized quark mass, when every term in the r.h.s. of the equation is computed at vanishing bare mass. We use the values for the critical mass provided in [35],

$$\begin{aligned} m_\mathrm{cr}^{(1)}&= -0.20255651209\,C_\mathrm{F}~~~(c_\mathrm{sw}=1),\nonumber \\ m_\mathrm{cr}^{(1)}&= -0.32571411742\,C_\mathrm{F}~~~(c_\mathrm{sw}=0), \end{aligned}$$

(3.29)

with $C_\mathrm{F}=(N^2-1)/(2N)$, and the (tree-level) value of the Sheikholeslami–Wohlert (SW) coefficient $c_\mathrm{sw}$ indicating whether the computation is performed with or without an $\text{ O }(a)$-improved action.

The entries of the renormalization matrix admit a similar expansion,

$$\begin{aligned} Z(g_0,L/a) = 1 + g_0^2 Z^{(1)}(L/a) + \mathcal{O}(g_0^4), \end{aligned}$$

(3.30)

where we have indicated explicitly the dependence of the quantities on the bare coupling and the lattice spacing-rescaled renormalization scale $a\mu =a/L$. The explicit expression of the one-loop order coefficient $Z^{(1)}$ for the multiplicatively renormalizable operators $\mathcal{Q}_1^\pm$ is

$$\begin{aligned} Z^{(1)}= & {} -\left\{ \frac{F^{(1)}}{F^{(0)}} + \frac{F^{(1)}_b}{F^{(0)}} + m_\mathrm{cr}^{(1)} \frac{\partial }{\partial m_0} \log F^{(0)} \right\} \nonumber \\&+\, \left( \frac{3}{2}-\alpha \right) \left[ \frac{f_1^{(1)}}{f_1^{(0)}} + \frac{f_{1;b}^{(1)}}{f_1^{(0)}} + m_\mathrm{cr}^{(1)}\frac{\partial }{\partial m_0}\log f_1^{(0)} \right] \nonumber \\&+\, \alpha \left[ \frac{k_1^{(1)}}{k_1^{(0)}} + \frac{k_{1;b}^{(1)}}{k_1^{(0)}} + m_\mathrm{cr}^{(1)}\frac{\partial }{\partial m_0} \log k_1^{(0)} \right] , \end{aligned}$$

(3.31)

while for the entries of each $2 \times 2$ submatrix that renormalizes operator pairs one has

$$\begin{aligned} Z_{ij}^{(1)}=-\mathcal{A}_{ik}^{(1)}\left[ \left( \mathcal{A}^{(0)} \right) ^{-1} \right] _{kj}, \end{aligned}$$

(3.32)

with

$$\begin{aligned} \mathcal{A}_{ij}^{(0)}=&\frac{F_{ij}^{(0)}}{\left[ f_1^{(0)}\right] ^{3/2-\alpha }\left[ k_1^{(0)} \right] ^{\alpha }},\nonumber \\ \mathcal{A}_{ij}^{(1)} =&\left\{ \left[ F_{ij}^{(1)} + F_{ij;b}^{(1)} + m_\mathrm{cr}^{(1)}\frac{\partial }{\partial m_0}F_{ij}^{(0)} \right] \right. \nonumber \\&-\left( \frac{3}{2}{-}\alpha \right) \left[ \frac{f_1^{(1)}}{f_1^{(0)}} {+} \frac{f_{1;b}^{(1)}}{f_1^{(0)}} + m_\mathrm{cr}^{(1)}\frac{\partial }{\partial m_0} \log f_1^{(0)} \right] F_{ij}^{(0)}\nonumber \\&-\alpha \left. \left[ \frac{k_1^{(1)}}{k_1^{(0)}} + \frac{k_{1;b}^{(1)}}{k_1^{(0)}} + m_\mathrm{cr}^{(1)}\frac{\partial }{\partial m_0} \log k_1^{(0)} \right] F_{ij}^{(0)} \right\} \nonumber \\&\times \left[ f_1^{(0)}\right] ^{\alpha -3/2}\left[ k_1^{(0)}\right] ^{-\alpha } . \end{aligned}$$

(3.33)

Contributions with the label “b” arise from the boundary terms that are needed in addition to the SW term in order to achieve full $\text{ O }(a)$ improvement of the action in the SF [33]. They obviously vanish in the unimproved case. We will set them to zero in the improved case as well, since they vanish in the continuum limit and thus will not contribute to our results for NLO anomalous dimensions.¹¹

Fig. 2

Feynman diagrams of the self-energy type entering the one-loop computation of $F_{k;s}^\pm$

Fig. 3

Feynman diagrams with gluon exchanges between quark lines entering the one-loop computation of $F_{k;s}^\pm$

The computation of the r.h.s. of the four-quark operator correlators $F_{k;s}^\pm$ requires the evaluation of the Feynman diagrams in Fig. 1 at tree level, and of those in Figs. 2 and 3 at one loop. The one-loop expansion of the boundary-to-boundary correlators $f_1$ and $k_1$ is meanwhile known from [36]. Each diagram can be written as a loop sum of a Dirac trace in time-momentum representation, where the Fourier transform is taken over space only. The sums have been performed numerically in double precision arithmetics using a Fortran 90 code, for all even lattice sizes ranging from $L/a=4$ to $L/a=48$. The results have been cross-checked against those of an independent C++ code, also employing double precision arithmetics.

The expected asymptotic expansion for the one-loop coefficient of renormalization constants is (operator and scheme indices not explicit)

$$\begin{aligned} Z^{(1)} (L/a) = \sum _{n=0}^\infty \left( \frac{a}{L}\right) ^n\left\{ r_n + s_n\,\ln (L/a) \right\} . \end{aligned}$$

(3.34)

In particular, the coefficient $s_0$ of the log that survives the continuum limit will be the corresponding entry of the anomalous dimension matrix, while the finite part $r_0$ will contribute to the one-loop matching coefficients we are interested in. In particular, one has

$$\begin{aligned}{}[\mathcal{X}_O^{(1)}]_\mathrm{SF;lat} = r_0, \end{aligned}$$

(3.35)

which is the required input for the matching condition in Eq. (3.10). We thus proceed as follows:

• (i) Compute tree-level and one-loop diagrams for all correlation functions.

• (ii) Construct one-loop renormalization constants using Eqs. (3.31) and (3.32).

• (iii) Fit the results to the ansatz in Eq. (3.34) as a function of (a / L), using the known value of the entries of the leading-order anomalous dimension matrix $\gamma ^{(0)}$ as fixed parameters, and extract $r_0$.

The description of the procedure employed to extract the finite parts as well as our results are provided in Appendix E.

3.5 NLO SF anomalous dimensions

Having collected $[\mathcal{X}_O^{(1)}]_\mathrm{SF;lat}$, $[\mathcal{X}_O^{(1)}]_\mathrm{cont;lat}$, $\gamma ^{(1);\mathrm cont}$ and $\mathcal{X}^{(1)}_\mathrm{g}$ we have finally been able to compute the matrix $\gamma ^{(1);\mathrm SF}$ for both the “+” and the “−” operator basis and for all the 18 schemes presented in Sect. 3.3. The results are collected in Appendix F.

We have performed two strong consistency checks of our calculation:

• In our one-loop perturbative computation, we have obtained $[\mathcal{X}_O^{(1)}]_\mathrm{SF;lat}$ for both $c_\mathrm{sw}=0$ and $c_\mathrm{sw}=1$ values. The results for $[\mathcal{X}_O^{(1)}]_\mathrm{cont;lat}$ are known for generic values of $c_\mathrm{sw}$. We have thus been able to compute $[\mathcal{X}_O^{(1)}]_\mathrm{SF;cont}$ for both $c_\mathrm{sw}=0$ and $c_\mathrm{sw}=1$ in such a way to check its independence from $c_\mathrm{sw}$.

• For the “+” operators, we have checked the independence of $\gamma ^{(1);\mathrm SF}$ from the reference scheme used (either the RI-MOM or the ${\overline{\mathrm{MS}}}$). This is a strong check of the calculations from the literature of the NLO anomalous dimensions $\gamma ^{(1);\mathrm cont}$ and one-loop matching coefficients $[\mathcal{X}_O^{(1)}]_\mathrm{cont;lat}$ in both the RI-MOM and ${\overline{\mathrm{MS}}}$ scheme.

The resulting values of $\gamma ^{(1)}$ exhibit a strong scheme dependence. In order to define a reference scheme for each operator, we have devised a criterion that singles out those schemes with the smallest NLO corrections: given the matrix

$$\begin{aligned} 16 \pi ^2\, \gamma ^{(1);\mathrm SF}\, (\gamma ^{(0)})^{-1}, \end{aligned}$$

(3.36)

we compute the trace and the determinant of each non-trivial submatrix, and look for the smallest absolute value of both quantities. Remarkably, in all cases (2–3 and 4–5 operator doublets, both in the Fierz $+$ and − sectors) this is satisfied by the scheme given by $s=6,~\alpha =3/2$.

4 Renormalization group running in perturbation theory

In this section we will discuss the perturbative computation of the RG running factor $\tilde{U}(\mu )$ in Eq. (2.40). The main purpose of this exercise is to understand the systematic of perturbative truncation, both in view of our own non-perturbative computation of the RG running factor [37] (which involves a matching to NLO perturbation theory around the electroweak scale), and in order to assess the extensive use of NLO RG running down to hadronic scales in the phenomenological literature. In view of our upcoming publication of a non-perturbative determination of the anomalous dimensions for QCD with $N_{\mathrm {\scriptstyle f}}=2$, the analysis below will be performed for that case; the qualitative conclusions are independent of the precise value of $N_{\mathrm {\scriptstyle f}}$. The scale will be fixed using the value $\Lambda _\mathrm{\scriptscriptstyle QCD}^{{\overline{\mathrm{MS}}};N_{\mathrm {\scriptstyle f}}=2}=310(20)~\mathrm{MeV}$ quoted in [38].

At leading order in perturbation theory the running factor is given by $U_\mathrm{\scriptscriptstyle LO}$ in Eq. (2.33). Beyond LO, the running factor is given by Eq. (2.34), where $W(\mu )$ satisfies Eq. (2.36). In the computation of W, the $\beta$ and the $\gamma$ functions are known only up to three loops and two loops, respectively. In order to asses the systematic, we will compute the running factor for several approximations that will be labeled through a pair of numbers “$n_\gamma /n_\beta$” where $n_\gamma$ is the order used for the $\gamma$ function while $n_\beta$ is the order used for the $\beta$ function. We will consider the following cases:

• (i) “1/1”, i.e. the LO approximation in which $W\equiv 1$;

• (ii) “2/2”, in which both $\gamma$ and $\beta$ are taken at NLO;

• (iii) “2/3”, in which $\beta$ is taken at NNLO and $\gamma$ at NLO;

• (iv) “+3/3”, in which $\beta$ is taken at NNLO and for the NNLO coefficient $\gamma _2$ we use a guesstimate given by $\gamma _2\gamma _1^{-1}=\gamma _1\gamma _0^{-1}$;

• (v) “−3/3”, in which $\beta$ is taken at the NNLO and for the NNLO coefficient $\gamma _2$ we use a guesstimate given by $\gamma _2\gamma _1^{-1}=-\gamma _1\gamma _0^{-1}$;

Beyond LO we have first computed the perturbative expansion of the running factor, Eqs. (2.34) and (2.41), by including all the $J_n$’s corresponding to the highest order used in the combinations of $\beta$/$\gamma$ functions chosen above. The $J_n$ have been computed from Eqs. (2.42) and (2.43) setting the unknown coefficients to zero. Namely: $J_1$ in the 2/2 case, $J_1$ and $J_2$ (with $\gamma _2=0$) in the 2/3 case, $J_1$ and $J_2$ with $\gamma _2$ set to the guesstimates above in the +3/3 and −3/3 cases. We have compared these results with the numerical solution of Eq. (2.36) in which the perturbative expansions for $\gamma$ and $\beta$ at the chosen orders are plugged in. We have chosen two cases in which perturbation theory seems particularly ill-behaved, namely the matrix for operators 4 and 5 with both Fierz + and − in the RI-MOM scheme, and we show the comparison in Fig. 4. As one can see, the two methods are not in very good agreement in the region of few GeV scales. This is obvious because, by expanding W in powers of $g^2$ and including only the first/second coefficients $J_1$, $J_2$, substantial information is lost.

We have then included in the perturbative expansion the next order, computed from Eqs. (2.43) and (2.44), setting again the unknown coefficients to zero. Namely: $J_2$ (with $b_2=\gamma _2=0$) in the 2/2 case, $J_3$ (with $b_3=\gamma _3=\gamma _2=0$) in the 2/3 case, $J_3$ (with $b_3=\gamma _3=0$ and $\gamma _2$ set to the guesstimates above) in the +3/3 and −3/3 cases. The comparison, again with the corresponding numerical solution of Eq. (2.36) (which remains unchanged), is shown in Fig. 5 and shows a reasonable agreement for the Fierz + matrix while there is still noticeable disagreement for some of the Fierz − matrix elements.

Fig. 4

RG running matrix for the Op 4, 5 in the RI scheme. Top half a Fierz $+$. Bottom half b Fierz −. The four cases $n_\gamma /n_\beta = \{2/2, 2/3, +3/3, -3/3\}$ are plotted, respectively, in red, black, magenta and blue. Dashed lines correspond to the numerical integration of $W(\mu )$. Solid lines correspond to the perturbative expansion up to ${\text{ O }}(g^2)$ (i.e. $J_1$) for the 2/2 case and up to ${\text{ O }}(g^4)$ (i.e. $J_2$) for the 2/3, $+3/3$ and $-3/3$ cases

Fig. 5

RG running matrix for the Op 4, 5 in the RI scheme. Top half a Fierz $+$. Bottom half b Fierz −. The four cases $n_\gamma /n_\beta = \{2/2, 2/3, +3/3, -3/3\}$ are plotted, respectively, in red, black, magenta and blue. Dashed lines correspond to the numerical integration of $W(\mu )$. Solid lines correspond to the perturbative expansion up to ${\text{ O }}(g^4)$ (i.e. $J_2$) for the 2/2 case and up to ${\text{ O }}(g^6)$ (i.e. $J_3$) for the 2/3, $+3/3$ and $-3/3$ cases

Fig. 6

RG running matrix for the Op 4, 5 Fierz − in the RI scheme. Top half a The four cases $n_\gamma /n_\beta = \{2/2, 2/3, +3/3, -3/3\}$ are plotted, respectively, in red, black, magenta and blue. Dashed lines correspond to the numerical integration of $W(\mu )$. Solid lines correspond to the perturbative expansion up to ${\text{ O }}(g^6)$ (i.e. $J_3$) for the 2/2 case and up to ${\text{ O }}(g^8)$ (i.e. $J_4$) for the 2/3, $+3/3$ and $-3/3$ cases. Bottom half b Comparison of the results for the numerical integration of $W(\mu )$ when matched at $\overline{g}^2_{{\overline{\mathrm{MS}}}}(M_P)$ with the perturbative expansion at the order used in Fig. 4 (solid lines) and Fig. 5 (dashed lines)

Table 1 Values for the RG running coefficients at $\mu =3\,\mathrm{GeV}$ for the four doublets of operators in three different schemes (${\overline{\mathrm{MS}}}$, RI and a chosen SF scheme). We quote here, as our best result, the case $n_\gamma /n_\beta$=2/3 obtained by numerical integration. The systematic errors have been estimated by computing the maximal deviation between the central value and the values of the 2/2, $+3/3$ and $-3/3$ numerical solutions

Table 2 Values for the RG running coefficients at $\mu =3\,\mathrm{GeV}$ for the four doublets of operators in three different schemes (${\overline{\mathrm{MS}}}$, RI and a chosen SF scheme). We quote here the 2/2 result from the perturbative expansion at $\text{ O }(g^2)$, which is the case usually considered in the literature, both for phenomenological application and in lattice computations. The systematic errors have been estimated by computing the maximal deviation between the central value and the values of the 2/2, 2/3, $+3/3$ and $-3/3$ numerical solutions. It is worth noticing the large asymmetric errors which occur in particular in the 45 Fierz $+$ and − matrices (especially in the RI scheme)

In the Fierz − case we have thus proceeded by introducing the next order, namely: $J_3$ (with $b_2=\gamma _2=b_3=\gamma _3=0$) in the 2/2 case, $J_4$ (with $b_4=b_3=\gamma _4=\gamma _3=\gamma _2=0$) in the 2/3 case, $J_4$ (with $b_4=b_3=\gamma _4=\gamma _3=0$ and $\gamma _2$ set to the guesstimates above) in the $+3/3$ and $-3/3$ cases. The comparison, again with the corresponding numerical solution of Eq. (2.36), is shown in Fig. 6a. The agreement between the numerical solution and the perturbative expansion further improves in all cases except for the 55 matrix element in the $\pm 3/3$ cases where the perturbative expansion further moves away from the numerical solution. From both examples of Fierz ± 4–5 matrix, we understand that by including more and more orders in the perturbative expansion of $W(\mu )$ Eq. (2.41), we approximate better and better the numerical solution of Eq. (2.36),¹² which can thus be considered the best approximation of the running factor given a fixed-order computation of the $\beta$ and $\gamma$ functions.

There is still a subtle technical issue concerning the numerical integration of Eq. (2.36) which needs to be discussed, because it becomes relevant in practice. Since $\gamma$ and $\beta$ have simple expressions in terms of $\overline{g}(\mu )$ rather than in terms of $\mu$, Eq. (2.36) is most easily solved by rewriting it in terms of the derivative with respect to the coupling, viz.

$$\begin{aligned} \widetilde{W}'(g) = -\widetilde{W}(g)\,\frac{\gamma (g)}{\beta (g)} +\frac{\gamma _0}{b_0g}\,\widetilde{W}(g), \end{aligned}$$

(4.1)

where $\widetilde{W}(\overline{g}(\mu )) \equiv W(\mu )$. While both terms on the right-hand side diverge as $g \rightarrow 0$, the divergence cancels in the sum due to Eq. (2.37). However, it is not straightforward to implement the latter initial condition at the level of the numerical solution to Eq. (4.1): a stable numerical solution requires fixing the initial condition Eq. (2.37) at an extremely small value of the coupling, and consequently the use of a sophisticated and computationally expensive integrator. A simpler solution is to substitute Eq. (2.37) by an initial condition of the form

$$\begin{aligned} \widetilde{W}(g_\mathrm{i}) = \mathbf {1} + g_\mathrm{i}^2 J_1 + g_\mathrm{i}^4 J_2 + g_\mathrm{i}^6 J_3 + g_\mathrm{i}^8 J_4 + \cdots , \end{aligned}$$

(4.2)

at some very perturbative coupling $g_\mathrm{i}$ (but still a significantly larger value than required by Eq. (2.37)), where we include exactly the same coefficients $J_n$, $n=1,\ldots$ that we use in the perturbative expansion of the running factor, and which are computed by using the same amount of perturbative information employed in the ratio $\gamma /\beta$ used for the numerical integration.¹³ Note that indeed the numerical value of $g_\mathrm{i}$ needs not be extremely small for this to make physical sense, e.g. for $N_{\mathrm {\scriptstyle f}}=2$ (which will be of particular interest to us) and at the Planck scale one has $\overline{g}^2_{{\overline{\mathrm{MS}}}}(M_P)\approx 0.221 \leftrightarrow \alpha _\mathrm{s}^{{\overline{\mathrm{MS}}}}(M_P) \approx 0.0176$ and $\overline{g}^2_\mathrm{SF}(M_P)$ differs with respect to $\overline{g}^2_{{\overline{\mathrm{MS}}}}(M_P)$ only on the third decimal digit.

In Fig. 6b we compare the results for the numerical integration of $W(\mu )$ when matched at $g_\mathrm{i}$ with the perturbative expansion at the order used in Figs. 4 and 5, respectively, and the results turn out to be indistinguishable. We have also changed the value of the coupling chosen for the matching in a broad range of $\overline{g}^2$ without observing any noticeable difference in the solution. These checks prove the stability of the numerical procedure and give us confidence in the corresponding results, which will be used below to assess the systematic uncertainties. In the following we will not consider anymore the perturbative expansion of the running factor except for the 2/2 case where only $J_1$ is included (we will call this 2 / 2 at $\text{ O }(g^2)$), which is the case usually considered in the literature, both for phenomenological application and in lattice computations.

According to the previous discussion, we have chosen to quote as our best estimate of the running factors the 2/3 results (which encode the maximum of information at our disposal for the $\beta$ and $\gamma$ functions) obtained through numerical integration. They are presented in Table 1 at the scale $\mu =3\,\mathrm{GeV}$. In alternative we quote also the results for the 2/2 case perturbatively expanded at $\text{ O }(g^2)$ (i.e. including $J_1$ only), which are the results usually considered in the literature. We present them in Table 2 again at the scale $\mu =3\,\mathrm{GeV}$.

The systematic uncertainties in Table 1 (respectively, Table 2) are estimated by considering the maximal deviation of the 2/3 case (respectively, the 2/2, $\text{ O }(g^2)$ case) from the other 3 (respectively, 4) numerical cases.

Fig. 7

RG running matrices for the Fierz $+$ Op. 2, 3 (top half) and Op. 4, 5 (bottom half) in the ${\overline{\mathrm{MS}}}$ scheme. Solid lines correspond to the LO plotted (cyan) and the perturbative expansion for the NLO 2/2 case up to ${\text{ O }}(g^2)$ – i.e. including $J_1$ (green). Dashed lines correspond to the numerical solution for $W(\mu )$ in the cases $n_\gamma /n_\beta = \{2/2, 2/3, +3/3, -3/3\}$, respectively, in red, black, magenta and blue

Fig. 8

RG running matrices for the Fierz $+$ Op. 2, 3 (top half) and Op. 4, 5 (bottom half) in the RI scheme. Solid lines correspond to the LO (cyan) and the perturbative expansion for the NLO 2/2 case up to ${\text{ O }}(g^2)$ – i.e. including $J_1$ (green). Dashed lines correspond to the numerical solution for $W(\mu )$ in the cases $n_\gamma /n_\beta = \{2/2, 2/3, +3/3, -3/3\}$, respectively, in red, black, magenta and blue

Fig. 9

RG running matrices for the Fierz $+$ Op. 2, 3 (top half) and Op. 4, 5 (bottom half) in the SF scheme. Solid lines correspond to the LO (cyan) and the perturbative expansion for the NLO 2/2 case up to ${\text{ O }}(g^2)$ – i.e. including $J_1$ (green). Dashed lines correspond to the numerical solution for $W(\mu )$ in the cases $n_\gamma /n_\beta = \{2/2, 2/3, +3/3, -3/3\}$, respectively, in red, black, magenta and blue

Fig. 10

RG running matrices for the Fierz − Op. 2, 3 (top half) and Op. 4, 5 (bottom half) in the ${\overline{\mathrm{MS}}}$ scheme. Solid lines correspond to the LO (cyan) and the perturbative expansion for the NLO 2/2 case up to ${\text{ O }}(g^2)$ – i.e. including $J_1$ (green). Dashed lines correspond to the numerical solution for $W(\mu )$ in the cases $n_\gamma /n_\beta = \{2/2, 2/3, +3/3, -3/3\}$, respectively, in red, black, magenta and blue

Fig. 11

RG running matrices for the Fierz − Op. 2, 3 (top half) and Op. 4, 5 (bottom half) in the RI scheme. Solid lines correspond to the LO (cyan) and the perturbative expansion for the NLO 2/2 case up to ${\text{ O }}(g^2)$ – i.e. including $J_1$ (green). Dashed lines correspond to the numerical solution for $W(\mu )$ in the cases $n_\gamma /n_\beta = \{2/2, 2/3, +3/3, -3/3\}$, respectively, in red, black, magenta and blue

Fig. 12

RG running matrices for the Fierz − Op. 2, 3 (top half) and Op. 4, 5 (bottom half) in the SF scheme. Solid lines correspond to the LO (cyan) and the perturbative expansion for the NLO 2/2 case up to ${\text{ O }}(g^2)$ – i.e. including $J_1$ (green). Dashed lines correspond to the numerical solution for $W(\mu )$ in the cases $n_\gamma /n_\beta = \{2/2, 2/3, +3/3, -3/3\}$, respectively, in red, black, magenta and blue

The results for the LO running factor $U_\mathrm{\scriptscriptstyle LO}(\mu )$ Eq. (2.33) and the numerically integrated $\tilde{U}(\mu )$ running factors beyond LO (ii)–(v) described above are illustrated in Figs. 7, 8, 9, 10, 11, and 12 together with the 2/2 $\text{ O }(g^2)$ perturbative expansion, for the four doublets of operators and three different schemes (${\overline{\mathrm{MS}}}$, RI and a chosen SF scheme).

Some important observations are:

• The convergence of LO respect to NLO and NNLO seems to be slow in all the schemes under investigation for almost all the operators. In particular, for the matrix elements involving tensor current (4–5 submatrices) the convergence is very poor. Note that the LO anomalous dimensions for these submatrices are already very large compared with the others.

• the 2/3 numerical running factors have always symmetric systematic errors, because most of the systematics is due to the inclusion of the guesstimate for $\gamma _2$ with + and − sign, and these effects turn out to be always symmetric with respect to the 2/3 (and also 2/2) cases.

• the 2/2 $\text{ O }(g^2)$ running factors are, for several matrix elements, quite far from the 2/3 (and also the 2/2) numerical ones. Possibly even further away than the $\pm 3/3$ and have thus very large, asymmetric errors.

• For both 4–5 submatrices (Fierz + and −) the ratio $\gamma _1\gamma _0^{-1}$ turns out to have large matrix elements. As a consequence, our plausibility argument for the guesstimates $\gamma _2\gamma _1^{-1}=\pm \gamma _1\gamma _0^{-1}$ leads to large systematic uncertainties. In particular, for the 54 matrix element the error is huge in the RI scheme and large also in the ${\overline{\mathrm{MS}}}$ and SF schemes, already for the 2/3 numerical solution (for the 2/2 $\text{ O }(g^2)$ perturbative expansion the situation is much worse). This obviously poses serious doubts on all the computations of $\Delta F=2$ matrix elements beyond the Standard Model which uses perturbative running (in all cases through the 2/2 $\text{ O }(g^2)$ expansion) down to scales of 3 $\mathrm{GeV}$ or less.

5 Conclusions

In this paper we have reviewed the renormalization and RG running properties of the four-quark operators relevant for BSM analyses, and introduced a family of SF schemes that allow one to compute them in a fully non-perturbative way. Our non-perturbative results for $N_{\mathrm {\scriptstyle f}}=2$ QCD will be presented in a separate publication [39].¹⁴ Here we have focused on the perturbative matching of our schemes to commonly used perturbative schemes and to RGI operators. One of our main results in this context is the full set of NLO operator anomalous dimensions in our SF schemes.

We have also conducted a detailed analysis of perturbative truncation effects in operator RG running in both SF schemes introduced here, and in commonly used ${\overline{\mathrm{MS}}}$ and RI-MOM schemes. We conclude that when NLO perturbation theory is used to run the operators from high-energy scales down to the few GeV range, large truncation effects appear. One striking example is the mixing of tensor–tensor and scalar–scalar operators, where all the available indications point to extremely large anomalous dimensions and very poor perturbative convergence. One important point worth stressing is that, in the computation of the running factor $W(\mu )$, the use of the truncated perturbative expansion in Eq. (2.41) leads to a significantly worse behavior than the numerical integration of Eq. (2.36) with the highest available orders for $\gamma$ and $\beta$.

A context where these findings might have an important impact is e.g. the computation of BSM contributions to neutral kaon mixing. At present, few computations of the relevant $\Delta S=2$ operators exist with dynamical fermions [41,42,43,44,45], all of which use perturbative RG running (and, in the case of [44], perturbative operator renormalization as well). There are substantial discrepancies between the various results in [41,42,43,44,45], which may be speculated to stem, at least in part, from perturbative truncation effects. Another possible contribution to the discrepancy is the delicate pole subtraction required in the RI-MOM scheme – indeed, results involving perturbative renormalization and non-perturbative renormalization constants in RI-SMOM schemes are consistent. At any rate, future efforts to settle this issue, as well as similar studies for $\Delta B=2$ amplitudes, should put a strong focus on non-perturbative renormalization.

Appendices

Constraints on anomalous dimensions from chiral symmetry

In section 5.3 of [1] the authors derive an identity between the renormalization matrices for $(\mathcal{Q}_2^+,\mathcal{Q}_3^+)$ and $(\mathcal{Q}_2^-,\mathcal{Q}_3^-)$, valid in the RI-MOM scheme considered in that paper. Here we discuss how such an identity can be derived from generic considerations based on chiral symmetry, and how (or, rather, under which conditions) it can be generalized to other renormalization schemes.

Let us consider a renormalized matrix element of the form $\langle f|\bar{Q}_k^\pm |i\rangle$, where $Q_k^\pm$ is a parity-even operator and $|i,f\rangle$ are stable hadron states with the same, well-defined parity. Simple examples would be the matrix elements of $\Delta F=2$ operators providing the hadronic contribution to $K^0$–$\bar{K}^0$ or $B^0$–$\bar{B}^0$ oscillation amplitudes (cf. Sect. 2). Bare matrix elements can be extracted from suitable three-point Euclidean correlation functions

$$\begin{aligned}&\langle \mathcal{O}_f(x)\,Q_k^\pm (0)\,\mathcal{O}_i(y)\rangle \nonumber \\&\quad = \frac{1}{\mathcal{Z}}\int D[\psi ]D[\bar{\psi }]D[A]\,e^{-S}\,\mathcal{O}_f(x)\,Q_k^\pm (0)\,\mathcal{O}_i(y)\nonumber \\ \end{aligned}$$

(A.1)

where $\mathcal{O}_{i,f}$ are interpolating operators for the external states $|i,f\rangle$. If we perform a change of fermion variables of the form

$$\begin{aligned} \psi \rightarrow \psi '=e^{i\gamma _5 T}\psi ,\quad \bar{\psi }\rightarrow \bar{\psi }'=\bar{\psi }e^{i\gamma _5 T}, \end{aligned}$$

(A.2)

where $\psi$ is a fermion field with $N_{\mathrm {\scriptstyle f}}$ flavor components and T is a traceless matrix acting on flavor space, this will induce a corresponding transformation $Q_k^\pm \rightarrow {Q_k'}^\pm$, $\mathcal{O}_{i,f}\rightarrow \mathcal{O}'_{i,f}$ of the involved composite operators. If the regularized theory employed to define the path integral preserves exactly the $\mathrm{SU}(N_{\mathrm {\scriptstyle f}})_\mathrm{A}$ axial chiral symmetry of the formal continuum theory, the equality $\langle \mathcal{O}_f(x)\,Q_k^\pm (0)\,\mathcal{O}_i(y)\rangle =\langle \mathcal{O}'_f(x) \,{Q_k'}^\pm (0)\,\mathcal{O}'_i(y)\rangle$ will hold exactly; otherwise, it will only hold upon renormalization and removal of the cutoff. At the level of matrix elements, one will then have

$$\begin{aligned} \langle f|\bar{Q}_k^\pm |i\rangle _{(\psi ,\bar{\psi })} = \langle f|\bar{Q_k'}^\pm |i\rangle _{(\psi ',\bar{\psi }')} , \end{aligned}$$

(A.3)

where the subscript remarks that the interpretation of the operator depends on the fermion variables used on each side of the equation. If the flavor matrix T is not traceless, the argument will still hold if the fermion fields entering composite operators are part of a valence sector, employed only for the purpose of defining suitable correlation functions.

The result in Eq. (A.2) is at the basis e.g. of the definition of twisted-mass QCD lattice regularizations, and is discussed in more detail in [2,3,4]. Indeed, the rotation in Eq. (A.2) will in general transform the mass term of the action. One crucial remark at this point is that, if a mass-independent renormalization scheme is used, renormalization constants for any given composite operator will be independent of which fermion variables are employed in the computation of the matrix element.

Let us now consider a particular case of Eq. (A.2) given by

$$\begin{aligned} T = \frac{\pi }{4}\left( \begin{array}{rrrr} 1 &{} 0 &{} 0 &{} 0 \\ 0 &{} 1 &{} 0 &{} 0 \\ 0 &{} 0 &{} 1 &{} 0 \\ 0 &{} 0 &{} 0 &{} -1 \end{array}\right) , \end{aligned}$$

(A.4)

where $\psi =(\psi _1,\psi _2,\psi _3,\psi _4)^T$ comprises the four, formally distinct flavors that enter $Q_k^\pm ,\mathcal{Q}_k^\pm$. Under this rotation, the ten operators of the basis in Eq. (2.1) transform as

$$\begin{aligned} Q_1^\pm\rightarrow & {} i\mathcal{Q}_1^\pm ,\nonumber \\ Q_2^\pm\rightarrow & {} -i\mathcal{Q}_2^\mp ,\nonumber \\ Q_3^\pm\rightarrow & {} i\mathcal{Q}_3^\mp ,\\ Q_4^\pm\rightarrow & {} i\mathcal{Q}_4^\pm ,\nonumber \\ Q_5^\pm\rightarrow & {} i\mathcal{Q}_5^\pm .\nonumber \end{aligned}$$

(A.5)

In the case of operators 1, 4, 5 the rotation is essentially trivial, in that it preserves Fierz ($2 \leftrightarrow 4$ exchange) eigenstates. However, in the rotation of operators 2, 3 the Fierz eigenvalue is exchanged. One thus has, at the level of renormalized matrix elements,

$$\begin{aligned} \langle f|\bar{\mathbf {Q}}^+(\mu )|i\rangle _{(\psi ,\bar{\psi })} = R\langle f|\bar{\varvec{{\mathcal{Q}}}}^-(\mu )|i\rangle _{(\psi ',\bar{\psi }')} , \end{aligned}$$

(A.6)

where $\mathbf {Q}^+=(Q_2^+,Q_3^+)^T$, ${\varvec{{\mathcal{Q}}}}^-=(\mathcal{Q}_2^-,\mathcal{Q}_3^-)^T$, and $R=-i\tau ^3$. In the latter expression we have written explicitly the renormalization scale $\mu$. If we now use the RG evolution operators discussed in Sect. 2 to run Eq. (A.6) to another scale $\mu '$, one then has (recall that the continuum anomalous dimensions of $Q_k^+$ and $\mathcal{Q}_k^+$ – respectively, $Q_k^-$ and $\mathcal{Q}_k^-$ – are the same)

$$\begin{aligned} \langle f|\bar{\mathbf {Q}}^+(\mu ')|i\rangle _{(\psi ,\bar{\psi })}&= U^+(\mu ',\mu )\langle f|\bar{\mathbf {Q}}^+(\mu )|i\rangle _{(\psi ,\bar{\psi })}\nonumber \\&= U^+(\mu ',\mu )R\langle f|\bar{\varvec{{\mathcal{Q}}}}^-(\mu )|i\rangle _{(\psi ',\bar{\psi }')}\nonumber \\&= U^+(\mu ',\mu )R [U^-(\mu ',\mu )]^{-1}\nonumber \\&\quad \times \langle f|\bar{\varvec{{\mathcal{Q}}}}^-(\mu ')|i\rangle _{(\psi ',\bar{\psi }')}\nonumber \\&= U^+(\mu ',\mu )R [U^-(\mu ',\mu )]^{-1}\nonumber \\&\quad \times R^{-1}\langle f|\bar{\mathbf {Q}}^+(\mu ')|i\rangle _{(\psi ,\bar{\psi })}, \end{aligned}$$

(A.7)

which implies

$$\begin{aligned} U^+(\mu ',\mu ) = R U^-(\mu ',\mu ) R^{-1}\quad \forall \mu ,\mu '. \end{aligned}$$

(A.8)

It is then immediate that the anomalous dimension matrices entering $U^\pm$ are related as

$$\begin{aligned} \gamma ^+ = \left( \begin{array}{cc} \gamma _{22}^+ &{} \gamma _{23}^+ \\ \gamma _{32}^+ &{} \gamma _{33}^+ \end{array}\right) = \tau ^3 \gamma ^- \tau ^3 = \left( \begin{array}{cc} \gamma _{22}^- &{} -\gamma _{23}^- \\ -\gamma _{32}^- &{} \gamma _{33}^- \end{array}\right) . \end{aligned}$$

(A.9)

The correct interpretation of this identity is that, given an anomalous dimension matrix for, say, $Q_{2,3}^+$ and $\mathcal{Q}_{2,3}^+$, one can use Eq. (A.9) to construct a correct anomalous dimension matrix for $Q_{2,3}^-$ and $\mathcal{Q}_{2,3}^-$, and vice versa. However, it does not guarantee that, given two different renormalization conditions for each fierzing, the resulting matrices of anomalous dimensions will satisfy Eq. (A.9). This will only be the case if the renormalization conditions can be related to each other by the rotation in Eq. (A.4); otherwise, the result of applying Eq. (A.9) to the $\gamma ^-$ that follows from the condition imposed on Fierz − operators will lead to value of $\gamma ^+$ in a different renormalization scheme than the one defined by the renormalization condition imposed directly on Fierz + operators.

The RI-MOM conditions of [1], as well as typical ${\overline{\mathrm{MS}}}$ renormalization conditions, result in schemes that satisfy the identity directly, since the quantities involved respect the underlying chiral symmetry – e.g. the amputated correlation functions used in RI-MOM rotate in a similar way to the three-point functions discussed above. Indeed, the known NLO anomalous dimensions in RI-MOM and ${\overline{\mathrm{MS}}}$ given in Appendix B, as well as (within uncertainties) the non-perturbative values of RI-MOM renormalization constants, fulfill Eq. (A.9). Our SF renormalization conditions, on the other hand, are not related among them via rotations with R, due to the chiral symmetry-breaking effects induced by the non-trivial boundary conditions imposed on the fields. As a consequence, the finite parts of the matrices of SF renormalization constants, and hence $\gamma _2^\mathrm{SF}$, do not satisfy the identity. It has to be stressed that, as a consequence of the existence of schemes where Eq. (A.9) is respected, the identity is satisfied by the universal matrices $\gamma _0^\pm$, as can be readily checked in Eq. (2.24); therefore, the violation of the identity in e.g. SF schemes appears only at $\mathcal{O}(g_0^4)$ in perturbation theory.

NLO anomalous dimensions in continuum schemes

The two-loop anomalous dimension matrices in the RI-MOM scheme (in Landau gauge) [5, 8] and ${\overline{\mathrm{MS}}}$ scheme [8] are given by (the factor $(4\pi )^{-4}$ has been omitted below to simplify the notation):

$$\begin{aligned} \gamma ^{+,(1);\mathrm{RI}}_{22}&= \frac{(297+16 \log (2)) N^2+45}{6 N^2}\\&\quad -N_f\frac{2 (15+4\log (2))}{3 N} ,\\ \gamma ^{+,(1);\mathrm{RI}}_{23}&= \frac{2 \left( 4 N^2 (45+2\log (2))-9\right) }{3 N}\\&\quad -N_f\frac{4}{3} (15+4\log (2)) ,\\ \gamma ^{+,(1);\mathrm{RI}}_{32}&= \frac{(53+160 \log (2)) N^2+108}{12 N}\\&\quad -N_f\frac{2}{3} (1+2\log (2)) , \end{aligned}$$

$$\begin{aligned} \gamma ^{+,(1);\mathrm{RI}}_{33}&= \frac{-379 N^4+5 (99+32 \log (2)) N^2+45}{6 N^2}\\&\quad +\,N_f\frac{2 \left( 13 N^2-4 \log (2)-15\right) }{3 N} ,\\ \end{aligned}$$

$$\begin{aligned} \gamma ^{+,(1);\mathrm{RI}}_{44}&= \frac{-379 N^4+2 (261-88 \log (2)) N^3+140 (3+2\log (2)) N^2-4 (-6+60 \log (2)) N-81}{6 N^2}\\&\quad +\,N_f\frac{2 \left( 13 N^2+(-15+8\log (2)) N-4\log (2)-15\right) }{3 N} ,\\ \gamma ^{+,(1);\mathrm{RI}}_{45}&= \frac{(157-368 \log (2)) N^3+(-494+556 \log (2)) N^2-4 (-39+30 \log (2)) N-72}{36 N^2}\\&\quad +\,N_f\frac{((-11+16 \log (2)) N-20 \log (2)+28)}{18 N} ,\\ \gamma ^{+,(1);\mathrm{RI}}_{54}&= \frac{4 \left( (-165+16 \log (2)) N^3+(-230+76 \log (2)) N^2-4 (-39+30 \log (2)) N-72\right) }{3 N^2}\\&\quad -\,N_f\frac{8 ((15+16 \log (2)) N-20 \log (2)+28)}{3 N} ,\\ \gamma ^{+,(1);\mathrm{RI}}_{55}&= \frac{343 N^4-2 (-343+616 \log (2)) N^3+4 (-95+142 \log (2)) N^2+(504+720 \log (2)) N-531}{18 N^2}\\&\quad -\,N_f\frac{2 \left( 13 N^2+(41-56 \log (2)) N+52 \log (2)-11\right) }{9 N} , \end{aligned}$$

$$\begin{aligned} \gamma ^{-,(1);\mathrm{RI}}_{22}&= \frac{15}{2 N^2}+\frac{8 \log (2)}{3}+\frac{99}{2}\\&\quad -\,N_f\frac{2 (15+4\log (2))}{3 N} ,\\ \gamma ^{-,(1);\mathrm{RI}}_{23}&= -\frac{8}{3} (45+2\log (2)) N\\&\quad +\,\frac{6}{N} +N_f\frac{4}{3} (15+4\log (2)) ,\\ \gamma ^{-,(1);\mathrm{RI}}_{32}&= -\frac{1}{12} (53+160 \log (2)) N\\&\quad -\,\frac{9}{N}+N_f\frac{2}{3} (1+2\log (2)) ,\\ \gamma ^{-,(1);\mathrm{RI}}_{33}&= \frac{-379 N^4+5 (99+32 \log (2)) N^2+45}{6 N^2}\\&\quad +\,N_f\frac{2 \left( 13 N^2-4 \log (2)-15\right) }{3 N}, \end{aligned}$$

$$\begin{aligned} \gamma ^{-,(1);\mathrm{RI}}_{44}&= \frac{-379 N^4+2 (-261+88 \log (2)) N^3+140 (3+2\log (2)) N^2+24 (-1+10\log (2)) N-81}{6 N^2}\\&\quad +\,N_f\frac{2 \left( 13 N^2-(-15+8\log (2)) N-4\log (2)-15\right) }{3 N} ,\\ \gamma ^{-,(1);\mathrm{RI}}_{45}&= \frac{(-157+368 \log (2)) N^3+(-494+556 \log (2)) N^2+12 (-13+10\log (2)) N-72}{36 N^2}\\&\quad -\,N_f\frac{((-11+16\log (2)) N+20\log (2)-28)}{18 N} ,\\ \gamma ^{-,(1);\mathrm{RI}}_{54}&= -\frac{4 \left( (-165+16 \log (2)) N^3+(230-76 \log (2)) N^2+(156-120 \log (2)) N+72\right) }{3 N^2}\\&\quad -\,N_f\frac{8 ((15+16 \log (2)) N+4 (-7+5\log (2)))}{3 N} ,\\ \gamma ^{-,(1);\mathrm{RI}}_{55}&= \frac{343 N^4+14 (-49+88 \log (2)) N^3+4 (-95+142 \log (2)) N^2-72 (7+10\log (2)) N-531}{18 N^2}\\&\quad -\,N_f\frac{2 \left( 13 N^2+(-41+56 \log (2)) N+52 \log (2)-11\right) }{9 N}, \end{aligned}$$

$$\begin{aligned} \gamma ^{+,(1);{{\overline{\mathrm{MS}}}}}_{22}&= \frac{15}{2 N^2}+\frac{137}{6}-N_f\frac{22}{3 N} ,\\ \gamma ^{+,(1);{{\overline{\mathrm{MS}}}}}_{23}&= \frac{200 N}{3}-\frac{6}{N}-N_f\frac{44}{3} ,\\ \gamma ^{+,(1);{{\overline{\mathrm{MS}}}}}_{32}&= \frac{71 N}{4}+\frac{9}{N}-N_f 2,\\ \gamma ^{+,(1);{{\overline{\mathrm{MS}}}}}_{33}&= -\frac{203 N^2}{6}+\frac{479}{6}+\frac{15}{2 N^2}+N_f\left( \frac{10 N}{3}-\frac{22}{3 N}\right) ,\\ \gamma ^{+,(1);{{\overline{\mathrm{MS}}}}}_{44}&= -\frac{203 N^2}{6}+\frac{107 N}{3}+\frac{136}{3}-\frac{12}{N}-\frac{107}{2 N^2}\\&\quad +\,N_f\left( \frac{10 N}{3}-\frac{2}{3}-\frac{10}{3 N}\right) ,\\ \gamma ^{+,(1);{{\overline{\mathrm{MS}}}}}_{45}&= -\frac{N}{36}-\frac{31}{9}+\frac{9}{N}-\frac{4}{N^2}+N_f\left( \frac{1}{9 N}-\frac{1}{18}\right) ,\\ \gamma ^{+,(1);{{\overline{\mathrm{MS}}}}}_{54}&= -\frac{364 N}{3}-\frac{704}{3}-\frac{208}{N}-\frac{320}{N^2}\\&\quad +\,N_f\left( \frac{136}{3}+\frac{176}{3 N}\right) ,\\ \gamma ^{+,(1);{{\overline{\mathrm{MS}}}}}_{55}&= \frac{343 N^2}{18}+21 N-\frac{188}{9}+\frac{44}{N}+\frac{21}{2 N^2}\\&\quad +\,N_f\left( -\frac{26 N}{9}-6+\frac{2}{9 N}\right) , \end{aligned}$$

$$\begin{aligned} \gamma ^{-,(1);{{\overline{\mathrm{MS}}}}}_{22}&= \frac{15}{2 N^2}+\frac{137}{6}-N_f\frac{22}{3 N} ,\\ \gamma ^{-,(1);{{\overline{\mathrm{MS}}}}}_{23}&= -\frac{200 N}{3}+\frac{6}{N}+N_f\frac{44}{3} ,\\ \gamma ^{-,(1);{{\overline{\mathrm{MS}}}}}_{32}&= -\frac{71 N}{4}-\frac{9}{N}+N_f 2 ,\\ \gamma ^{-,(1);{{\overline{\mathrm{MS}}}}}_{33}&= -\frac{203 N^2}{6}+\frac{479}{6}+\frac{15}{2 N^2}+N_f\left( \frac{10 N}{3}-\frac{22}{3 N}\right) , \end{aligned}$$

$$\begin{aligned}\gamma ^{-,(1);{{\overline{\mathrm{MS}}}}}_{44}&= -\frac{203 N^2}{6}-\frac{107 N}{3}+\frac{136}{3}+\frac{12}{N}-\frac{107}{2 N^2}\\&\quad +\,N_f\left( \frac{10 N}{3}+\frac{2}{3}-\frac{10}{3 N}\right) ,\\ \gamma ^{-,(1);{{\overline{\mathrm{MS}}}}}_{45}&= \frac{N}{36}-\frac{31}{9}-\frac{9}{N}-\frac{4}{N^2}+N_f \left( \frac{1}{18}+\frac{1}{9 N}\right) ,\\ \gamma ^{-,(1);{{\overline{\mathrm{MS}}}}}_{54}&= \frac{364 N}{3}-\frac{704}{3}+\frac{208}{N}-\frac{320}{N^2}\\&\quad +\,N_f\left( \frac{176}{3 N}-\frac{136}{3}\right) ,\\ \gamma ^{-,(1);{{\overline{\mathrm{MS}}}}}_{55}&= \frac{343 N^2}{18}-21 N-\frac{188}{9}-\frac{44}{N}+\frac{21}{2 N^2}\\&\quad +\,N_f\left( -\frac{26 N}{9}+6+\frac{2}{9 N}\right) . \end{aligned}$$

Perturbative expansion of RG evolution for $N_{\mathrm {\scriptstyle f}}=3$

It is well known that the condition in Eq. (2.42) that determines the leading non-trivial coefficient in the NLO perturbative expansion of the RG evolution operator, Eq. (2.41), is ill-behaved for the operators $Q_{2,3}^\pm$ for $N_{\mathrm {\scriptstyle f}}=30$ and, more relevantly, for $N_{\mathrm {\scriptstyle f}}=3$ [46, 47]. The reason is that, when Eq. (2.42) is written as a linear system, the $4 \times 4$ matrix that multiplies the vector of elements of $J_1$ has zero determinant, rendering the system indeterminate.

A simple way to understand the anatomy of this problem in greater detail proceeds by writing the explicit solution to Eq. (2.42) as a function of the parameter $\epsilon =3-N_{\mathrm {\scriptstyle f}}$; if the NLO anomalous dimension matrix in the scheme under consideration is written as

$$\begin{aligned} \gamma _1^\pm = \frac{1}{(4\pi )^4}\left( \begin{array}{cc} g_{22}^\pm &{} g_{23}^\pm \\ g_{32}^\pm &{} g_{33}^\pm \\ \end{array} \right) \end{aligned}$$

(C.1)

then one finds

$$\begin{aligned} J_1^\pm = \frac{1}{\epsilon }\,J_{1,-1}^\pm + J_{1,0}^\pm + \epsilon J_{1,1}^\pm + \mathcal{O}(\epsilon ^2), \end{aligned}$$

(C.2)

with

$$\begin{aligned}&J_{1,-1}^\pm = \frac{1}{(4\pi )^2}\left( \begin{array}{cc} 0 &{} \pm \frac{1}{2}(g_{22}^\pm - g_{33}^\pm ) - \frac{3}{4}g_{23}^\pm + \frac{1}{3}g_{32}^\pm \\ 0 &{} 0 \end{array}\right) ,\end{aligned}$$

(C.3)

$$\begin{aligned}&J_{1,0}^\pm = \frac{1}{(4\pi )^2}\nonumber \\&\quad \times \,\left( \begin{array}{cc} \frac{1}{162}(128-9g_{22}^\pm \mp 3g_{32}^\pm ) &{} \frac{1}{27}(\pm 128 \mp g_{22}^\pm -g_{32}^\pm \pm g_{33}^\pm )\\ -\frac{1}{36}g_{32}^\pm &{} \frac{1}{162}(-1024 \pm 3g_{32}^\pm -9g_{33}^\pm ) \\ \end{array}\right) ,\nonumber \\ \end{aligned}$$

(C.4)

$$\begin{aligned}&J_{1,1}^\pm = \frac{1}{(4\pi )^2}\nonumber \\&\quad \times \,\left( \begin{array}{cc} \frac{1}{4374}(172+18g_{22}^\pm \pm 9g_{32}^\pm ) &{} \pm \frac{1}{2187}( 516+6g_{22}^\pm \pm 7g_{32}^\pm -6g_{33}^\pm )\\ \frac{1}{972}g_{32}^\pm &{} \frac{1}{4374}(-1376 \mp 9g_{32}^\pm +18g_{33}^\pm )\\ \end{array}\right) . \end{aligned}$$

(C.5)

In the limit $\epsilon \rightarrow 0$ the element 23 of $J_1^\pm$ diverges; it is easy to check that the aforementioned $4 \times 4$ matrix, consistently, has determinant $\propto \epsilon$. A similar expansion of the matrices $\tilde{U}_\mathrm{\scriptscriptstyle LO}^\pm (\mu ) \equiv [\alpha _{\mathrm{\scriptscriptstyle s}}(\mu )]^{-\frac{\gamma _0^\pm }{2b_0}}$ yields

$$\begin{aligned} \tilde{U}_\mathrm{\scriptscriptstyle LO}^\pm (\mu ) = \tilde{U}_\mathrm{\scriptscriptstyle LO,0}^\pm (\mu ) + \epsilon \tilde{U}_\mathrm{\scriptscriptstyle LO,1}^\pm (\mu ) + \mathcal{O}(\epsilon ^2), \end{aligned}$$

(C.6)

with

$$\begin{aligned} \tilde{U}_\mathrm{\scriptscriptstyle LO,0}^\pm (\mu )&= \left( \begin{array}{cc} \alpha _{\mathrm{\scriptscriptstyle s}}^{-1/9}(\mu ) &{} \mp \frac{2}{3}[\alpha _{\mathrm{\scriptscriptstyle s}}^{8/9}(\mu )-\alpha _{\mathrm{\scriptscriptstyle s}}^{-1/9}(\mu )] \\ 0 &{} \alpha _{\mathrm{\scriptscriptstyle s}}^{8/9}(\mu )\\ \end{array}\right) , \end{aligned}$$

(C.7)

$$\begin{aligned} \tilde{U}_\mathrm{\scriptscriptstyle LO,1}^\pm (\mu )&= \left( \begin{array}{cc} \frac{2}{243}\alpha _{\mathrm{\scriptscriptstyle s}}^{-1/9}(\mu ) &{} \pm \frac{4}{729}[8\alpha _{\mathrm{\scriptscriptstyle s}}^{8/9}(\mu )+\alpha _{\mathrm{\scriptscriptstyle s}}^{-1/9}(\mu )] \\ 0 &{} -\frac{16}{243}\alpha _{\mathrm{\scriptscriptstyle s}}^{8/9}(\mu )\\ \end{array}\right) \nonumber \\&\quad \times \log [\alpha _{\mathrm{\scriptscriptstyle s}}(\mu )]. \end{aligned}$$

(C.8)

When these expressions are plugged in Eq. (2.34), and the perturbative expansion Eq. (2.41) is used, one gets

$$\begin{aligned} \tilde{U}^\pm (\mu )= & {} ~ \tilde{U}_\mathrm{\scriptscriptstyle LO,0}^\pm (\mu ) + \overline{g}^2(\mu )\left[ \frac{1}{\epsilon } \tilde{U}_\mathrm{\scriptscriptstyle LO,0}^\pm (\mu ) J_{1,-1}^\pm \right. \nonumber \\&\left. +\, \tilde{U}_\mathrm{\scriptscriptstyle LO,0}^\pm (\mu ) J_{1,0}^\pm + \tilde{U}_\mathrm{\scriptscriptstyle LO,1}^\pm (\mu ) J_{1,-1}^\pm + \mathcal{O}(\epsilon ) \right] \nonumber \\&+\, \mathcal{O}(\overline{g}^4(\mu )), \end{aligned}$$

(C.9)

which is still divergent as $\epsilon \rightarrow 0$. This implies, in particular, that RGI operators cannot be defined consistently using the above form of the perturbative expansion for W. The RG evolution operator $U(\mu _2,\mu _1)=[\tilde{U}(\mu _2)]^{-1} \tilde{U}(\mu _1)$, on the other hand, is finite: the divergent part has the form

$$\begin{aligned}&\frac{1}{\epsilon }\left\{ \overline{g}^2(\mu _1)U_\mathrm{\scriptscriptstyle LO}^\pm (\mu _2,\mu _1) J_{1,-1}^\pm - \overline{g}^2(\mu _2) J_{1,-1}^\pm U_\mathrm{\scriptscriptstyle LO}^\pm (\mu _2,\mu _1) \right\} \nonumber \\&\quad = \frac{\pm \frac{1}{2}(g_{22}^\pm - g_{33}^\pm ) - \frac{3}{4}g_{23}^\pm + \frac{1}{3}g_{32}^\pm }{4\pi \epsilon } \,M, \end{aligned}$$

(C.10)

with

$$\begin{aligned} M= & {} U_\mathrm{\scriptscriptstyle LO}^\pm (\mu _2,\mu _1) \left( \begin{array}{cc} 0 &{} \alpha _{\mathrm{\scriptscriptstyle s}}(\mu _2) \\ 0 &{} 0 \end{array}\right) \nonumber \\&- \left( \begin{array}{cc} 0 &{} \alpha _{\mathrm{\scriptscriptstyle s}}(\mu _1) \\ 0 &{} 0 \end{array}\right) U_\mathrm{\scriptscriptstyle LO}^\pm (\mu _2,\mu _1) , \end{aligned}$$

(C.11)

and it is easy to check, using the explicit expression for $\tilde{U}_\mathrm{\scriptscriptstyle LO,0}^\pm (\mu )$ and the identity $U_\mathrm{\scriptscriptstyle LO}^\pm (\mu _2,\mu _1)=[U_\mathrm{\scriptscriptstyle LO,0}^\pm (\mu _2)]^{-1}U_\mathrm{\scriptscriptstyle LO,0}^\pm (\mu _1)$, that $M=0$.¹⁵ The full expression for $U(\mu _2,\mu _1)$ in the $\epsilon \rightarrow 0$ limit still receives contributions from $J_{1,-1}$, via the products with the $\mathcal{O}(\epsilon )$ terms in the expansion of $\tilde{U}_\mathrm{\scriptscriptstyle LO}$, which actually give rise to the only dependence of the expanded $U(\mu _2,\mu _1)$ on $\gamma _1^\pm$.

A number of solutions to this problem have been proposed in the literature [46,47,48,49], consisting of various regularization schemes to treat the singular terms in $3-N_{\mathrm {\scriptstyle f}}$. Here we note that the problem can be entirely bypassed by using the numerical integration of the RG equation in Eq. (2.36), as done in this paper to explore the case $N_{\mathrm {\scriptstyle f}}=2$ in detail. Indeed, applying exactly the same procedure for $N_{\mathrm {\scriptstyle f}}=3$ – i.e., solving Eq. (2.36) after having substituted the perturbative expressions for $\gamma$ and $\beta$ to any prescribed order – is well behaved numerically, which in turn allows one to construct both the RG evolution matrix and the RGI operators without trouble. The only point in the procedure where the expansion coefficient $J_1$ may enter explicitly is the initial condition in Eq. (4.2), where for the $N_{\mathrm {\scriptstyle f}}=2$ case we have employed $W(\mu _0)=\mathbf {1}+\overline{g}^2(\mu _0)J_1$ at some very high energy scale $\mu _0$. However, this can be replaced by the initial condition $W(\mu _0)=\mathbf {1}$ at an even higher scale, thus again avoiding the appearance of any singularity; it turns out that the required value of $\overline{g}^2(\mu )$ has to be extremely small, such that the systematics associated to the choice of coupling for the initial condition is negligible at the level of the result run down to $\overline{g}^2(\mu ) \sim 2$. This in turn requires using an expensive numerical integrator to work across several orders of magnitude, which is easy e.g. using standard Mathematica functions, provided proper care is taken to choose a stable integrator.

Fig. 13

W as a function of the coupling constant, with $N_{\mathrm {\scriptstyle f}}=3$ for operators $O_2$, $O_3$ fierz “$+$”, in the ${\overline{\mathrm{MS}}}$ scheme

As a crosscheck of the robustness of our numerical approach we have computed explicitly the function $W(\mu )$ for $N_{\mathrm {\scriptstyle f}}=3$, using our numerical integration and $W=\mathbf {1}$ as an initial condition, set at an extremely small value of the coupling. Our result for W, displayed in Fig. 13, can then be fitted to an ansatz where J is taken to have a polynomial dependence in $\overline{g}^2$, to check whether the first coefficient $J_1$ is compatible within systematic fit errors (obtained by trying different polynomial orders up to $\mathcal {O}(\overline{g}^8)$ and coupling values for the initial condition) with the one quoted in Eq. (2.30) of [48]. Note that in order to have a direct comparison it has to be taken into account that we are using a different relative normalization between operators $O_2$, $O_3$ than the one in [48] which translates into a factor $-2$ and $-1/2$, respectively, for $[J_1^\mathrm{fit}]_{23}$, $[J_1^\mathrm{fit}]_{32}$ and that, since we are working with the renormalization constants instead of the Wilson coefficients, the convention used for the J is this work corresponds to $J^T$ in [48]. What we obtain is

$$\begin{aligned} J_1^\mathrm{fit} = \frac{1}{(4\pi )^2}\left( \begin{array}{cc} -1.0470(8) &{} 70.13(38)\\ -1.39583(1) &{} 5.78550(8) \end{array} \right) , \end{aligned}$$

(C.12)

which is indeed well compatible with the above-mentioned result. Note that the coefficient 23 contains some precise numerical value of the parameter t employed in [48] to regularize the divergence of J in $3-N_{\mathrm {\scriptstyle f}}$.

As a further crosscheck, we have also compared the result of computing the $N_{\mathrm {\scriptstyle f}}=2$ evolution with the two possible initial conditions. The outcome is that, if the value of the coupling at which $W=\mathbf {1}$ is sufficiently small, the two results are equal up to several significant figures down to values of the coupling $\overline{g}^2 \gtrsim 2$, where the hadronic regime is entered.

Finite parts of RI-MOM renormalization constants in Landau gauge

In this appendix we gather the results for the finite part of the one-loop matching coefficients $[\mathcal{X}_O^{(1)}]_\mathrm{RI;lat}$ between the lattice and the RI-MOM scheme in Landau gauge. They can be extracted from [11] and are given by

$$\begin{aligned}{}[\mathcal{X}_{22}^{+,(1)}]_\mathrm{RI;lat}= & {} 0.0272369 \,c_\mathrm{sw}^2+0.0485167 \,c_\mathrm{sw}-0.294894 \nonumber ,\\ [\mathcal{X}_{23}^{+,(1)}]_\mathrm{RI;lat}= & {} 0.0218485 \,c_\mathrm{sw}^2+0.0632421 \,c_\mathrm{sw}+0.0753979 \nonumber ,\\ [\mathcal{X}_{32}^{+,(1)}]_\mathrm{RI;lat}= & {} 0.00755569 \nonumber ,\\ [\mathcal{X}_{33}^{+,(1)}]_\mathrm{RI;lat}= & {} -0.00553581 \,c_\mathrm{sw}^2-0.0463464 \,c_\mathrm{sw}-0.362656 ,\nonumber \\ [\mathcal{X}_{44}^{+,(1)}]_\mathrm{RI;lat}= & {} 0.00538842 \,c_\mathrm{sw}^2-0.0147254 \,c_\mathrm{sw}-0.351294 \nonumber ,\\ [\mathcal{X}_{45}^{+,(1)}]_\mathrm{RI;lat}= & {} 0.000303451 \,c_\mathrm{sw}^2+0.000878362 \,c_\mathrm{sw}-0.00178318 \nonumber ,\\ [\mathcal{X}_{54}^{+,(1)}]_\mathrm{RI;lat}= & {} -0.0728282 \,c_\mathrm{sw}^2-0.210807 \,c_\mathrm{sw}-0.266293 \nonumber ,\\ [\mathcal{X}_{55}^{+,(1)}]_\mathrm{RI;lat}= & {} 0.0442301 \,c_\mathrm{sw}^2+0.0977049 \,c_\mathrm{sw}-0.290267,\nonumber \\ \end{aligned}$$

(D.1)

$$\begin{aligned}{}[\mathcal{X}_{22}^{-,(1)}]_\mathrm{RI;lat}= & {} 0.0272369 \,c_\mathrm{sw}^2+0.0485167 \,c_\mathrm{sw}-0.294894 \nonumber ,\\ [\mathcal{X}_{23}^{-,(1)}]_\mathrm{RI;lat}= & {} -0.0218485 \,c_\mathrm{sw}^2-0.0632421 \,c_\mathrm{sw}-0.0753979\nonumber ,\\ [\mathcal{X}_{32}^{-,(1)}]_\mathrm{RI;lat}= & {} -0.00755569\nonumber ,\\ [\mathcal{X}_{33}^{-,(1)}]_\mathrm{RI;lat}= & {} -0.00553581 \,c_\mathrm{sw}^2-0.0463464 \,c_\mathrm{sw}-0.362656,\nonumber \\ [\mathcal{X}_{44}^{-,(1)}]_\mathrm{RI;lat}= & {} -0.01646 \,c_\mathrm{sw}^2-0.0779674 \,c_\mathrm{sw}-0.374019\nonumber ,\\ [\mathcal{X}_{45}^{-,(1)}]_\mathrm{RI;lat}= & {} -0.00151725 \,c_\mathrm{sw}^2-0.00439181 \,c_\mathrm{sw}+0.0013602\nonumber ,\\ [\mathcal{X}_{54}^{-,(1)}]_\mathrm{RI;lat}= & {} 0.0145656 \,c_\mathrm{sw}^2+0.0421614 \,c_\mathrm{sw}+0.24599\nonumber ,\\ [\mathcal{X}_{55}^{-,(1)}]_\mathrm{RI;lat}= & {} 0.0223817 \,c_\mathrm{sw}^2+0.0344629 \,c_\mathrm{sw}-0.257729.\nonumber \\ \end{aligned}$$

(D.2)

Table 3 Numerical results of the one-loop finite parts for operators $\text {VA--AV, PS--SP}$ in the 18 SF renormalization schemes under investigation defined by the source s and the parameter $\alpha$ as in Eq. (3.33)

Table 4 Numerical results of the one-loop finite parts for operators $\text {VA--AV, PS--SP}$ in the 18 SF renormalization schemes under investigation defined by the source s and the parameter $\alpha$ as in Eq. (3.33). These results have been computed including the clover term in the fermionic action

Table 5 Numerical results of the one-loop finite parts for operators $\text {PS + SP, T}{\tilde{T}}$ in the 18 SF renormalization schemes under investigation defined by the source s and the parameter $\alpha$ as in Eq. (3.33)

Table 6 Numerical results of the one-loop finite parts for operators $\text {PS + SP, T}{\tilde{T}}$ in the 18 SF renormalization schemes under investigation defined by the source s and the parameter $\alpha$ as in Eq. (3.33). These results have been computed including the clover term in the fermionic action

Finite parts of SF renormalization constants

In this appendix we discuss how to determine the dependence on a / L of the one-loop renormalization constants $Z^{(1)}$ defined in Sect. 3. The approach is essentially an application of the present context to the techniques discussed in Appendix D of [50].

Defining $\ell =L/a$ we will hence consider $F(\ell )=Z^{(1)}$ as a pure function of $\ell =\{\ell _1,\dots ,\ell _N\}$. We will also assume that all divergences have been removed from $F(\ell )$, which in general means linear divergences related to the additive renormalization of quark masses and proportional to the one-loop value of the critical mass $m_{cr}^{(1)}$, and logarithmic divergences proportional to a LO anomalous dimension. To ensure the robustness of our method we performed separate fits, and we checked, for each ansatz, the fitted value of $\gamma ^{(0)}$ was the correct one within the available precision. First of all a roundoff error has to be assigned to $F(\ell )$, which takes into account the uncertainties coming from the numerical computation itself. Following [50], we choose as an estimate for this error, in the case that the computation has been carried out in double precision,

$$\begin{aligned} \delta F(\ell ) \equiv \epsilon (\ell ) |F(\ell )| , \end{aligned}$$

(E.1)

$$\begin{aligned} \epsilon (\ell )=\left( \frac{\ell }{2} \right) ^3 \times 10^{-14}. \end{aligned}$$

(E.2)

As showed in Sect. 3 the expected behavior of $F(\ell )$ leads to the consideration of an asymptotic expansion of the form

$$\begin{aligned} F(\ell )=\alpha _0+\sum _{k=1}^{n}\frac{1}{\ell ^k}(\alpha _k+\beta _k\log \ell )+R_n(\ell ) , \end{aligned}$$

(E.3)

where the residue $R_n(\ell )$ is expected to decrease faster as $\ell \rightarrow \infty$ than any of the terms in the sum. To determine the coefficients $(\alpha _k,\beta _k)$ we minimize a quadratic form in the residues

$$\begin{aligned} \chi ^2=(F-f\xi )^T(F-f\xi ), \end{aligned}$$

(E.4)

where F and $\xi$ are the $N-$ and $(2n+1)-$column vectors $(F(\ell _1),\dots ,F(\ell ))^T$ and $(\alpha _0,\alpha _1,\dots ,\alpha _n,\beta _1,\dots ,\beta _n)^T$, respectively, and f is the $N\times (2n+1)$ matrix

$$\begin{aligned} f=\begin{pmatrix} 1 &{} \ell _1^{-1} &{} \cdots &{} \ell _1^{-n} &{} \ell _1^{-1}\log \ell _1 &{} \cdots &{} \ell _1^{-n}\log \ell _1\\ 1 &{} \ell _2^{-1} &{} \cdots &{} \ell _2^{-n} &{} \ell _2^{-1}\log \ell _2 &{} \cdots &{} \ell _2^{-n}\log \ell _2\\ \vdots &{} \vdots &{} \vdots &{} \vdots &{} \vdots &{} \vdots &{} \vdots \\ 1 &{} \ell _N^{-1} &{} \cdots &{} \ell _N^{-n} &{} \ell _N^{-1}\log \ell _N &{} \cdots &{} \ell _N^{-n}\log \ell _N \end{pmatrix} \end{aligned}$$

(E.5)

Again following [50], we have not introduced a matrix of weights in the definition of $\chi ^2$. A necessary condition to minimize $\chi ^2$ is

$$\begin{aligned} f\xi =PF \end{aligned}$$

(E.6)

where we have assumed that the columns of f are linearly independent vectors (assuming $2n+1 \ll N$), and P is the projector onto the subspace of $\mathbb {R}^N$ generated by them. Equation (E.6) can be solved using the singular value decomposition of f, which has the form of

$$\begin{aligned} f=USV^T \end{aligned}$$

(E.7)

where U is an $N \times (2n+1)$ matrix such that

$$\begin{aligned} U^TU=\mathbf {1},\quad UU^T=P \end{aligned}$$

(E.8)

S is diagonal, and $(2n+1)\times (2n+1)$ matrix V is orthonormal. With this decomposition one has

$$\begin{aligned} \xi =VS^{-1}U^TF \, . \end{aligned}$$

(E.9)

Finally, the uncertainty in the result for $\xi _k$ can be modeled using error propagation as

$$\begin{aligned} \delta \xi _k^2=\sum _{l=1}^N(VS^{-1}U^T)_{kl}^2(\delta F)_l^2 , \end{aligned}$$

(E.10)

where $(\delta F)_k \equiv \delta F(\ell _k)$.

As a remark on the above method regarding practical applications, it has to be pointed out that the choice of Eq. (E.4) for the quadratic form $\chi ^2$ implies, in particular, that small values of $\ell$ might be given excessive weight. This problem has been dealt with by considering a range $[\ell _\mathrm{min},\ell _\mathrm{max}]$ with changing $\ell _\mathrm{min}$. For this work the better convergence in results for $(\alpha _k,\beta _k)$ was given by $\ell _\mathrm{min}=16$ and $\ell _\mathrm{max}=46$. The estimation of systematic uncertainty of the fitting procedure has be performed using the proposal by the authors of [50]. We considered two independent fits at order n and $n+1$, i.e. extending the Ansatz in Eq. (E.3) by terms $1/\ell ^{n+1}$ and $\log \ell / \ell ^{n+1}$ with coefficients $\alpha _{n+1}$ and $\beta _{n+1}$, respectively. The systematic uncertainty of the finite part $r_0=\alpha _0$ is defined as the difference of the value of the parameter $\alpha _0$ extracted by the two different fits. In the present work we have used $n=2$ in the fit Ansatz for the $\mathcal {O}$(a)-improved data, and $n=3$ for unimproved ones (Tables 3, 4, 5, 6).

NLO anomalous dimension matrices in SF schemes

$\alpha$	s	${\gamma ^{(1)}}_{23}^{+}$
0	1	$\begin{pmatrix}0.001519(10) + N_f[-0.000057850(2)] &{} 0.00983(2) + N_f[-0.00034710(1)] \\ 0.006188(1) + N_f[-0.000080203] &{} -0.006776(8) + N_f[-0.00001842(2)]\end{pmatrix}$
	2	$\begin{pmatrix}0.001080(8) + N_f[-0.000057850(2)] &{} 0.00936(2) + N_f[-0.00034710(1)] \\ 0.005504(3) + N_f[-0.000080203] &{} -0.006855(7) + N_f[-0.00001842(2)]\end{pmatrix}$
	3	$\begin{pmatrix}-0.001673(4) + N_f[-0.000057850(2)] &{} 0.00870(2) + N_f[-0.00034710(1)] \\ 0.008552(4) + N_f[-0.000080203] &{} -0.00651(1) + N_f[-0.00001842(2)]\end{pmatrix}$
	4	$\begin{pmatrix}0.000936(10) + N_f[-0.000057850(2)] &{} 0.00743(2) + N_f[-0.00034710(1)] \\ 0.007320(1) + N_f[-0.000080203] &{} -0.00290(2) + N_f[-0.00001842(2)]\end{pmatrix}$
	5	$\begin{pmatrix}0.000497(6) + N_f[-0.000057850(2)] &{} 0.00695(2) + N_f[-0.00034710(1)] \\ 0.0066351(1) + N_f[-0.000080203] &{} -0.00297(2) + N_f[-0.00001842(2)]\end{pmatrix}$
	6	$\begin{pmatrix}-0.002256(4) + N_f[-0.000057850(2)] &{} 0.00629(2) + N_f[-0.00034710(1)] \\ 0.009684(7) + N_f[-0.000080203] &{} -0.00263(2) + N_f[-0.00001842(2)]\end{pmatrix}$
3/2	1	$\begin{pmatrix}-0.000022(3) + N_f[-0.000057850(2)] &{} 0.00983(2) + N_f[-0.00034710(1)] \\ 0.006188(1) + N_f[-0.000080203] &{} -0.00832(2) + N_f[-0.00001842(2)]\end{pmatrix}$
	2	$\begin{pmatrix}-0.000461(1) + N_f[-0.000057850(2)] &{} 0.00936(2) + N_f[-0.00034710(1)] \\ 0.005504(3) + N_f[-0.000080203] &{} -0.00840(2) + N_f[-0.00001842(2)]\end{pmatrix}$
	3	$\begin{pmatrix}-0.003214(6) + N_f[-0.000057850(2)] &{} 0.00870(3) + N_f[-0.00034710(1)] \\ 0.008552(4) + N_f[-0.000080203] &{} -0.00805(3) + N_f[-0.00001842(2)]\end{pmatrix}$
	4	$\begin{pmatrix}-0.000605(3) + N_f[-0.000057850(2)] &{} 0.00743(1) + N_f[-0.00034710(1)] \\ 0.007320(1) + N_f[-0.000080203] &{} -0.004438(2) + N_f[-0.00001842(2)]\end{pmatrix}$
	5	$\begin{pmatrix}-0.0010440(6) + N_f[-0.000057850(2)] &{} 0.00695(1) + N_f[-0.00034710(1)] \\ 0.0066351(1) + N_f[-0.000080203] &{} -0.004516(2) + N_f[-0.00001842(2)]\end{pmatrix}$
	6	$\begin{pmatrix}-0.003797(7) + N_f[-0.000057850(2)] &{} 0.00629(2) + N_f[-0.00034710(1)] \\ 0.009684(7) + N_f[-0.000080203] &{} -0.004167(4) + N_f[-0.00001842(2)]\end{pmatrix}$
1	1	$\begin{pmatrix}0.000492(5) + N_f[-0.000057850(2)] &{} 0.00983(2) + N_f[-0.00034710(1)] \\ 0.006188(1) + N_f[-0.000080203] &{} -0.00780(2) + N_f[-0.00001842(2)]\end{pmatrix}$
	2	$\begin{pmatrix}0.000053(3) + N_f[-0.000057850(2)] &{} 0.00936(2) + N_f[-0.00034710(1)] \\ 0.005504(3) + N_f[-0.000080203] &{} -0.00788(2) + N_f[-0.00001842(2)]\end{pmatrix}$
	3	$\begin{pmatrix}-0.002700(4) + N_f[-0.000057850(2)] &{} 0.00870(2) + N_f[-0.00034710(1)] \\ 0.008552(4) + N_f[-0.000080203] &{} -0.00753(2) + N_f[-0.00001842(2)]\end{pmatrix}$
	4	$\begin{pmatrix}-0.000092(5) + N_f[-0.000057850(2)] &{} 0.00743(1) + N_f[-0.00034710(1)] \\ 0.007320(1) + N_f[-0.000080203] &{} -0.003924(7) + N_f[-0.00001842(2)]\end{pmatrix}$
	5	$\begin{pmatrix}-0.000530(2) + N_f[-0.000057850(2)] &{} 0.00695(1) + N_f[-0.00034710(1)] \\ 0.0066351(1) + N_f[-0.000080203] &{} -0.004002(7) + N_f[-0.00001842(2)]\end{pmatrix}$
	6	$\begin{pmatrix}-0.003283(5) + N_f[-0.000057850(2)] &{} 0.00629(2) + N_f[-0.00034710(1)] \\ 0.009684(7) + N_f[-0.000080203] &{} -0.003654(6) + N_f[-0.00001842(2)]\end{pmatrix}$

$\alpha$	s	${\gamma ^{(1)}}_{23}^{-}$
0	1	$\begin{pmatrix}0.00051(1) + N_f[-0.000118203(2)] &{} -0.00803(4) + N_f[0.00070922(1)] \\ -0.006546(2) + N_f[0.000063796] &{} -0.00660(1) + N_f[0.00056285(2)]\end{pmatrix}$
	2	$\begin{pmatrix}-0.002017(7) + N_f[-0.000118203(2)] &{} -0.00577(2) + N_f[0.00070922(1)] \\ -0.008464(2) + N_f[0.000063796] &{} -0.004564(4) + N_f[0.00056285(2)]\end{pmatrix}$
	3	$\begin{pmatrix}-0.002036(5) + N_f[-0.000118203(2)] &{} -0.00609(2) + N_f[0.00070922(1)] \\ -0.0069803(3) + N_f[0.000063796] &{} -0.006138(8) + N_f[0.00056285(2)]\end{pmatrix}$
	4	$\begin{pmatrix}-0.00049(1) + N_f[-0.000118203(2)] &{} -0.00498(2) + N_f[0.00070922(1)] \\ -0.007802(1) + N_f[0.000063796] &{} -0.00229(2) + N_f[0.00056285(2)]\end{pmatrix}$
	5	$\begin{pmatrix}-0.003025(6) + N_f[-0.000118203(2)] &{} -0.00272(1) + N_f[0.00070922(1)] \\ -0.009719(5) + N_f[0.000063796] &{} -0.00026(3) + N_f[0.00056285(2)]\end{pmatrix}$
	6	$\begin{pmatrix}-0.003045(4) + N_f[-0.000118203(2)] &{} -0.00305(1) + N_f[0.00070922(1)] \\ -0.008236(3) + N_f[0.000063796] &{} -0.00183(2) + N_f[0.00056285(2)]\end{pmatrix}$
3/2	1	$\begin{pmatrix}-0.001026(6) + N_f[-0.000118203(2)] &{} -0.00803(4) + N_f[0.00070922(1)] \\ -0.006546(2) + N_f[0.000063796] &{} -0.00814(3) + N_f[0.00056285(2)]\end{pmatrix}$
	2	$\begin{pmatrix}-0.003558(2) + N_f[-0.000118203(2)] &{} -0.00577(2) + N_f[0.00070922(1)] \\ -0.008464(2) + N_f[0.000063796] &{} -0.00611(2) + N_f[0.00056285(2)]\end{pmatrix}$
	3	$\begin{pmatrix}-0.003577(2) + N_f[-0.000118203(2)] &{} -0.00609(2) + N_f[0.00070922(1)] \\ -0.0069803(3) + N_f[0.000063796] &{} -0.00768(2) + N_f[0.00056285(2)]\end{pmatrix}$
	4	$\begin{pmatrix}-0.002034(3) + N_f[-0.000118203(2)] &{} -0.00498(1) + N_f[0.00070922(1)] \\ -0.007802(1) + N_f[0.000063796] &{} -0.003835(2) + N_f[0.00056285(2)]\end{pmatrix}$
	5	$\begin{pmatrix}-0.004566(5) + N_f[-0.000118203(2)] &{} -0.002722(8) + N_f[0.00070922(1)] \\ -0.009719(5) + N_f[0.000063796] &{} -0.00180(1) + N_f[0.00056285(2)]\end{pmatrix}$
	6	$\begin{pmatrix}-0.004586(5) + N_f[-0.000118203(2)] &{} -0.00305(1) + N_f[0.00070922(1)] \\ -0.008236(3) + N_f[0.000063796] &{} -0.003374(6) + N_f[0.00056285(2)]\end{pmatrix}$
*1	1	$\begin{pmatrix}-0.000512(8) + N_f[-0.000118203(2)] &{} -0.00803(4) + N_f[0.00070922(1)] \\ -0.006546(2) + N_f[0.000063796] &{} -0.00763(2) + N_f[0.00056285(2)]\end{pmatrix}$
	2	$\begin{pmatrix}-0.003044(2) + N_f[-0.000118203(2)] &{} -0.00577(2) + N_f[0.00070922(1)] \\ -0.008464(2) + N_f[0.000063796] &{} -0.00559(1) + N_f[0.00056285(2)]\end{pmatrix}$
	3	$\begin{pmatrix}-0.0030635(3) + N_f[-0.000118203(2)] &{} -0.00609(2) + N_f[0.00070922(1)] \\ -0.0069803(3) + N_f[0.000063796] &{} -0.00717(2) + N_f[0.00056285(2)]\end{pmatrix}$
	4	$\begin{pmatrix}-0.001521(6) + N_f[-0.000118203(2)] &{} -0.00498(1) + N_f[0.00070922(1)] \\ -0.007802(1) + N_f[0.000063796] &{} -0.003321(7) + N_f[0.00056285(2)]\end{pmatrix}$
	5	$\begin{pmatrix}-0.004052(3) + N_f[-0.000118203(2)] &{} -0.002722(8) + N_f[0.00070922(1)] \\ -0.009719(5) + N_f[0.000063796] &{} -0.00129(2) + N_f[0.00056285(2)]\end{pmatrix}$
	6	$\begin{pmatrix}-0.004072(3) + N_f[-0.000118203(2)] &{} -0.00305(1) + N_f[0.00070922(1)] \\ -0.008236(3) + N_f[0.000063796] &{} -0.00286(1) + N_f[0.00056285(2)]\end{pmatrix}$

$\alpha$	s	${\gamma ^{(1)}}_{45}^{+}$
0	1	$\begin{pmatrix}0.002303(3) + N_f[0.00012884(1)] &{} -0.00169802(7) + N_f[0.0000026054(2)] \\ 0.00172(8) + N_f[0.00179861(4)] &{} 0.00081(2) + N_f[-0.00035752(1)]\end{pmatrix}$
	2	$\begin{pmatrix}0.002685(2) + N_f[0.00012884(1)] &{} -0.0018770(3) + N_f[0.0000026054(2)] \\ 0.00336(7) + N_f[0.00179861(4)] &{} -0.00001(1) + N_f[-0.00035752(1)]\end{pmatrix}$
	3	$\begin{pmatrix}0.002077(4) + N_f[0.00012884(1)] &{} -0.0015930(3) + N_f[0.0000026054(2)] \\ 0.00232(7) + N_f[0.00179861(4)] &{} 0.00046(1) + N_f[-0.00035752(1)]\end{pmatrix}$
	4	$\begin{pmatrix}0.00558(2) + N_f[0.00012884(1)] &{} -0.0019027(6) + N_f[0.0000026054(2)] \\ 0.00795(6) + N_f[0.00179861(4)] &{} 0.00040(1) + N_f[-0.00035752(1)]\end{pmatrix}$
	5	$\begin{pmatrix}0.00597(2) + N_f[0.00012884(1)] &{} -0.002082(1) + N_f[0.0000026054(2)] \\ 0.00959(4) + N_f[0.00179861(4)] &{} -0.00043(1) + N_f[-0.00035752(1)]\end{pmatrix}$
	6	$\begin{pmatrix}0.00536(1) + N_f[0.00012884(1)] &{} -0.0017978(6) + N_f[0.0000026054(2)] \\ 0.00856(4) + N_f[0.00179861(4)] &{} 0.000048(10) + N_f[-0.00035752(1)]\end{pmatrix}$
3/2	1	$\begin{pmatrix}0.00076(2) + N_f[0.00012884(1)] &{} -0.00169802(7) + N_f[0.0000026054(2)] \\ 0.00172(8) + N_f[0.00179861(4)] &{} -0.000727(3) + N_f[-0.00035752(1)]\end{pmatrix}$
	2	$\begin{pmatrix}0.00114(1) + N_f[0.00012884(1)] &{} -0.0018770(3) + N_f[0.0000026054(2)] \\ 0.00336(7) + N_f[0.00179861(4)] &{} -0.001556(3) + N_f[-0.00035752(1)]\end{pmatrix}$
	3	$\begin{pmatrix}0.00054(2) + N_f[0.00012884(1)] &{} -0.0015930(3) + N_f[0.0000026054(2)] \\ 0.00232(7) + N_f[0.00179861(4)] &{} -0.001082(3) + N_f[-0.00035752(1)]\end{pmatrix}$
	4	$\begin{pmatrix}0.004044(3) + N_f[0.00012884(1)] &{} -0.0019027(5) + N_f[0.0000026054(2)] \\ 0.00795(3) + N_f[0.00179861(4)] &{} -0.001138(2) + N_f[-0.00035752(1)]\end{pmatrix}$
	5	$\begin{pmatrix}0.004426(5) + N_f[0.00012884(1)] &{} -0.002082(1) + N_f[0.0000026054(2)] \\ 0.00959(2) + N_f[0.00179861(4)] &{} -0.001966(5) + N_f[-0.00035752(1)]\end{pmatrix}$
	6	$\begin{pmatrix}0.003817(3) + N_f[0.00012884(1)] &{} -0.0017978(5) + N_f[0.0000026054(2)] \\ 0.00856(2) + N_f[0.00179861(4)] &{} -0.001492(5) + N_f[-0.00035752(1)]\end{pmatrix}$
1	1	$\begin{pmatrix}0.00128(1) + N_f[0.00012884(1)] &{} -0.00169802(7) + N_f[0.0000026054(2)] \\ 0.00172(8) + N_f[0.00179861(4)] &{} -0.000214(7) + N_f[-0.00035752(1)]\end{pmatrix}$
	2	$\begin{pmatrix}0.00166(1) + N_f[0.00012884(1)] &{} -0.0018770(3) + N_f[0.0000026054(2)] \\ 0.00336(7) + N_f[0.00179861(4)] &{} -0.001042(3) + N_f[-0.00035752(1)]\end{pmatrix}$
	3	$\begin{pmatrix}0.00105(1) + N_f[0.00012884(1)] &{} -0.0015930(3) + N_f[0.0000026054(2)] \\ 0.00232(7) + N_f[0.00179861(4)] &{} -0.000568(3) + N_f[-0.00035752(1)]\end{pmatrix}$
	4	$\begin{pmatrix}0.004557(7) + N_f[0.00012884(1)] &{} -0.0019027(6) + N_f[0.0000026054(2)] \\ 0.00795(4) + N_f[0.00179861(4)] &{} -0.000624(6) + N_f[-0.00035752(1)]\end{pmatrix}$
	5	$\begin{pmatrix}0.004940(9) + N_f[0.00012884(1)] &{} -0.002082(1) + N_f[0.0000026054(2)] \\ 0.00959(3) + N_f[0.00179861(4)] &{} -0.001453(2) + N_f[-0.00035752(1)]\end{pmatrix}$
	6	$\begin{pmatrix}0.004331(7) + N_f[0.00012884(1)] &{} -0.0017978(5) + N_f[0.0000026054(2)] \\ 0.00856(2) + N_f[0.00179861(4)] &{} -0.000979(1) + N_f[-0.00035752(1)]\end{pmatrix}$

$\alpha$	s	${\gamma ^{(1)}}_{45}^{-}$
0	1	$\begin{pmatrix}-0.01620(2) + N_f[0.00069509(2)] &{} 0.001678(1) + N_f[-0.0000064430(8)] \\ 0.00879(7) + N_f[-0.002125401(8)] &{} 0.00907(1) + N_f[-0.0004158697(7)]\end{pmatrix}$
	2	$\begin{pmatrix}-0.01676(2) + N_f[0.00069509(2)] &{} 0.001156(2) + N_f[-0.0000064430(8)] \\ 0.00887(6) + N_f[-0.002125401(8)] &{} 0.008779(9) + N_f[-0.0004158697(7)]\end{pmatrix}$
	3	$\begin{pmatrix}-0.01560(1) + N_f[0.00069509(2)] &{} 0.002266(2) + N_f[-0.0000064430(8)] \\ 0.00841(6) + N_f[-0.002125401(8)] &{} 0.008299(9) + N_f[-0.0004158697(7)]\end{pmatrix}$
	4	$\begin{pmatrix}-0.01236(2) + N_f[0.00069509(2)] &{} 0.001914(1) + N_f[-0.0000064430(8)] \\ 0.00755(4) + N_f[-0.002125401(8)] &{} 0.00896(1) + N_f[-0.0004158697(7)]\end{pmatrix}$
	5	$\begin{pmatrix}-0.01292(2) + N_f[0.00069509(2)] &{} 0.0013917(10) + N_f[-0.0000064430(8)] \\ 0.00762(3) + N_f[-0.002125401(8)] &{} 0.008664(8) + N_f[-0.0004158697(7)]\end{pmatrix}$
	6	$\begin{pmatrix}-0.01176(2) + N_f[0.00069509(2)] &{} 0.002502(3) + N_f[-0.0000064430(8)] \\ 0.00716(3) + N_f[-0.002125401(8)] &{} 0.008184(8) + N_f[-0.0004158697(7)]\end{pmatrix}$
3/2	1	$\begin{pmatrix}-0.01774(3) + N_f[0.00069509(2)] &{} 0.001678(1) + N_f[-0.0000064430(8)] \\ 0.00879(7) + N_f[-0.002125401(8)] &{} 0.007530(5) + N_f[-0.0004158697(7)]\end{pmatrix}$
	2	$\begin{pmatrix}-0.01830(4) + N_f[0.00069509(2)] &{} 0.001156(2) + N_f[-0.0000064430(8)] \\ 0.00887(6) + N_f[-0.002125401(8)] &{} 0.007238(3) + N_f[-0.0004158697(7)]\end{pmatrix}$
	3	$\begin{pmatrix}-0.01714(3) + N_f[0.00069509(2)] &{} 0.002266(2) + N_f[-0.0000064430(8)] \\ 0.00841(6) + N_f[-0.002125401(8)] &{} 0.006758(3) + N_f[-0.0004158697(7)]\end{pmatrix}$
	4	$\begin{pmatrix}-0.013901(2) + N_f[0.00069509(2)] &{} 0.0019137(8) + N_f[-0.0000064430(8)] \\ 0.00755(3) + N_f[-0.002125401(8)] &{} 0.007415(4) + N_f[-0.0004158697(7)]\end{pmatrix}$
	5	$\begin{pmatrix}-0.014461(4) + N_f[0.00069509(2)] &{} 0.0013917(6) + N_f[-0.0000064430(8)] \\ 0.00762(2) + N_f[-0.002125401(8)] &{} 0.007123(2) + N_f[-0.0004158697(7)]\end{pmatrix}$
	6	$\begin{pmatrix}-0.013305(5) + N_f[0.00069509(2)] &{} 0.002502(2) + N_f[-0.0000064430(8)] \\ 0.00716(2) + N_f[-0.002125401(8)] &{} 0.006643(2) + N_f[-0.0004158697(7)]\end{pmatrix}$
1	1	$\begin{pmatrix}-0.01722(3) + N_f[0.00069509(2)] &{} 0.001678(1) + N_f[-0.0000064430(8)] \\ 0.00879(7) + N_f[-0.002125401(8)] &{} 0.008043(7) + N_f[-0.0004158697(7)]\end{pmatrix}$
	2	$\begin{pmatrix}-0.01778(3) + N_f[0.00069509(2)] &{} 0.001156(2) + N_f[-0.0000064430(8)] \\ 0.00887(6) + N_f[-0.002125401(8)] &{} 0.007751(5) + N_f[-0.0004158697(7)]\end{pmatrix}$
	3	$\begin{pmatrix}-0.01663(2) + N_f[0.00069509(2)] &{} 0.002266(2) + N_f[-0.0000064430(8)] \\ 0.00841(6) + N_f[-0.002125401(8)] &{} 0.007271(5) + N_f[-0.0004158697(7)]\end{pmatrix}$
	4	$\begin{pmatrix}-0.013388(7) + N_f[0.00069509(2)] &{} 0.0019137(10) + N_f[-0.0000064430(8)] \\ 0.00755(3) + N_f[-0.002125401(8)] &{} 0.007928(6) + N_f[-0.0004158697(7)]\end{pmatrix}$
	5	$\begin{pmatrix}-0.013947(4) + N_f[0.00069509(2)] &{} 0.0013917(6) + N_f[-0.0000064430(8)] \\ 0.00762(3) + N_f[-0.002125401(8)] &{} 0.007637(4) + N_f[-0.0004158697(7)]\end{pmatrix}$
	6	$\begin{pmatrix}-0.01279(1) + N_f[0.00069509(2)] &{} 0.002502(2) + N_f[-0.0000064430(8)] \\ 0.00716(3) + N_f[-0.002125401(8)] &{} 0.007156(4) + N_f[-0.0004158697(7)]\end{pmatrix}$