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

We discuss the renormalization properties of the full set of ΔF=2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta F=2$$\end{document} 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 ΔF=1,2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta F=1,2$$\end{document} 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 MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\overline{\mathrm{MS}}}$$\end{document} and RI-MOM schemes. Large truncation effects are found for some of the operators considered.


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 O( 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 nonperturbative computation of the low-energy renormalization and RG running of all possible four-quark operators with net flavor change, by introducing appropriate Schrödinger a e-mail: david.preti@csic.es Functional (SF) renormalization schemes. In order to connect them with standard 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 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 F = 2 four-quark operators, stressing how considering it also allows one to obtain the anomalous dimensions of 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 fig-ures are collected after the main text, and a many technical details are discussed in appendices.

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 2 1 , and the labeling is adopted V → γ μ , A → γ μ γ 5 , S → 1, P → γ 5 , T → σ μν , T → 1 2 ε μνρτ σ ρτ , with σ μν = i 2 [γ μ , γ ν ]. In the above expression round parentheses indicate spin and color scalars, and subscripts are flavor labels. Note that operators Q ± k are parityeven, and Q ± k are parity-odd. It is important to stress that this framework is fairly general. For instance, with the assignments the operators Q − k vanish, while Q + 1 enters the SM amplitude for K 0 -K 0 mixing, and Q + 2,...,5 the contributions to the same amplitude from arbitrary extensions of the SM. Idem for B 0 (s) -B 0 (s) mixing with ψ 1 = ψ 3 = b, ψ 2 = ψ 4 = d/s. (2.4) If one instead chooses the assignments ψ 1 = s, ψ 2 = d, ψ 3 = ψ 4 = u, c, the resulting Q ± 1 will be the operators in the SM 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 and similarly in the parity-odd sector. If chiral symmetry is preserved by the regularization, both and D vanish. The main result of [1] is that D = 0 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 Q ± VA+AV that mediate Standard Model-allowed transitions renormalize multiplicatively, while operators that appear as effective interactions in extensions of the Standard Model do always mix. 1 Interestingly, in [1] some identities are derived that relate the renormalization matrices for (Q + 2 , Q + 3 ) and (Q − 2 , 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.

Callan-Symanzik equations
Theory parameters and operators are renormalized at the renormalization scale μ. The scale dependence of renormalized quantities is then governed by renormalization group evolution. We will consider QCD with N f quark flavors and N colors. The Callan-Symanzik equations satisfied by the gauge coupling and quark masses are of the form q ∂ ∂q g(q) = β(g(q)), (2.8) q ∂ ∂q m f (q) = τ (g(q))m f (q), (2.9) respectively, and satisfy the initial conditions 10) where f is a flavor label. Mass independence of the scheme is reflected in the fact that the beta function and mass anomalous dimension τ depend on the coupling and the number of flavors, but not on quark masses. Asymptotic perturbative expansions read (2.13) The universal coefficients of the perturbative beta function and mass anomalous dimension are 14) N .
We will deal with Euclidean correlation functions of gauge-invariant composite operators. Without loss of generality, let us consider correlation functions of the form which, expanding the total derivative, leads to where γ is a matrix of anomalous dimensions describing the mixing of {O k }, andγ l is the anomalous dimension of O l . For completeness, we have included a term which takes into account the dependence on the gauge parameter λ in covariant gauges; this term is absent in schemes like 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 β λ is given by and its perturbative expansion has the form (2.19) where the universal coefficient is given by In the Landau gauge (λ = 0) the term with β λ 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 β λ in all equations. From now on, in order to simplify the notation we will use the shorthand notation 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 22) and the perturbative expansion of the anomalous dimension matrix γ as 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 γ +,(0) 11

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 (μ 2 , μ 1 ) that evolves renormalized operators between the scales 3 μ 1 and μ 2 < μ 1 , By substituting into Eq. (2.21) one has the equation for U (μ 2 , μ 1 ) . the matrix product on the r.h.s.) with initial condition U (μ 1 , μ 1 ) = 1. Following a standard procedure, this differential equation for U can be converted into a Volterra-type integral equation and solved iteratively, viz.
where as usual the notation 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 3 Restricting the evolution operator to run towards the IR avoids unessential algebraic technicalities below. The running towards the UV can be trivially obtained by taking [U (μ 2 , μ 1 )] −1 .
T exp yielding the identity 4 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 μ 2 → ∞ ⇔ g(μ 2 ) → 0, and therefore so is the l.h.s. Thus, we define the RGI operator insertion as upon which we have an explicit expression to retrieve the RGI operator from the renormalized one at any value of the renormalization scale μ, provided the anomalous dimension and the beta function are known for scales ≥ μ, (2.32) Starting from the latter equation, it is easy to check explicitly thatÔ is invariant under a change of renormalization scheme. Note that the crucial step in the manipulation has been to add and subtract the term γ 0 b 0 g 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 γ (g)/β(g) ∼ 1/g. When γ is a matrix of anomalous dimensions this step becomes non-trivial, since in general [γ (g), γ 0 ] = 0; the derivation has thus to be performed somewhat more carefully.

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 γ (g)/β(g) the T-exponential becomes a standard exponential, since [γ 0 g 2 1 , γ 0 g 2 2 ] = 0 ∀g 1 , g 2 . One can then perform the integral trivially and write (2.33) When next-to-leading order corrections are included the Texponential becomes non-trivial. In order to make contact with the literature (see e.g. [5,8]), we write 5 (2.35) 5 The property underlying this equation is that the evolution operator can actually be factorized, in full generality, as (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 LO ; in particular, Note also that in the absence of mixing Eq. (2.36) can be solved explicitly to get (using W (0) = 1) Now it is easy, by analogy with the non-mixing case, to define RGI operators. We rewrite Eq. (2.25) as where 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 μ 1 → ∞ (resp. μ 2 → ∞), which implies that the vector of RGI operators can be obtained: W (μ) ≈ 1 + g 2 (μ)J 1 + g 4 (μ)J 2 + g 6 (μ)J 3 + g 8 (μ)J 4 + · · · (2.41) we find for the first four orders in the expansion the conditions (2.45) Modulo sign and normalization conventions (involving powers of 4π related to expanding in g 2 rather than α/(4π)), 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 γ . 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 g 2 in Eq. (2.41) up to zeroth, first, second, and third order, respectively.

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 The scheme-change factors X can be expanded perturbatively as 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, . (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 γ (0) . The relation between nextto-leading order coefficients for quark masses and operator anomalous dimensions are given by 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.

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 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 γ (1) is known. γ (1);MS can be found in [8][9][10], while γ (1);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 MS or RI prescription is applied to the latticeregulated theory. One can then proceed in two steps: (i) Compute the matching coefficient [X (1) O ] 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. MS or RI), [X (1) O ] cont;lat , and obtain the matching coefficient that enters Eq. (3.9) as (3.10) The one-loop matching coefficients [X (1) O ] 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 MS scheme they can be extracted from [12][13][14]). 7 We gather the RI-MOM results in Landau gauge in Appendix D. χ (1) g can be found in [15].

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 fourquark operators, the first applications involved the multiplica- 7 We are grateful to S. Sharpe for having converted for us, in the case of Fierz + operators, the MS scheme used in [12] to the one defined in [8]. tively renormalizable operators Q ± 1 of Eq. (2.1) (which, as explained above, enter Standard Model effective Hamiltonians for F = 1 and F = 2 processes) [28][29][30][31], as well as generic 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, where α, β are flavor indices, unprimed (primed) fields live at the x 0 = 0 (x 0 = T ) boundary, and is a spin matrix that must anticommute with γ 0 , so that the boundary fermion field does not vanish. This is a consequence of the structure of the conditions imposed on boundary fields, 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 parityodd 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 14) where S s is one of the five source operators and similarly for S s . The constant η k is a sign that ensures 8 We will also use the two-point functions of boundary sources, (3.23) Finally, we define the ratios where α is an arbitrary real parameter. The structure of F ± k;s and f 1 , k 1 is illustrated in Fig. 1.
We then proceed to impose renormalization conditions at bare coupling g 0 and scale μ = 1/L by generalizing the condition introduced in [28,29] for the renormalizable multiplicative operators Q ± 1 : the latter reads , (3.25) while, for operators that mix in doublets, we impose 9 , (3.26) and similarly for Q ± 4,5 . 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 , and therefore Z ± k;s,α 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 θ = 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 α. The latter can in principle take any value, but we will restrict our study to the choices α = 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 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 = s 2 . In practice, out of the ten possible schemes one is only left with six, viz. 10 It has to be stressed that this property is independent of the choice of θ and α. Thus, we are left with a total of 15 schemes for Q ± 1 , and 18 for each of the pairs (Q ± 2 , Q ± 3 ) and (Q ± 4 , Q ± 5 ).

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 where X is one of F ± k;s (x 0 ), f 1 , k 1 , or some combination thereof; where m 0 is the bare quark mass; and m (1) cr the oneloop 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], , and the (tree-level) value of the Sheikholeslami-Wohlert (SW) coefficient c sw indicating whether the computation is performed with or without an O(a)-improved action. The entries of the renormalization matrix admit a similar expansion, where we have indicated explicitly the dependence of the quantities on the bare coupling and the lattice spacingrescaled renormalization scale aμ = a/L. The explicit expression of the one-loop order coefficient Z (1) for the multiplicatively renormalizable operators Q ± 1 is while for the entries of each 2×2 submatrix that renormalizes operator pairs one has Contributions with the label "b" arise from the boundary terms that are needed in addition to the SW term in order to achieve full 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. 11 The computation of the r.h.s. of the four-quark operator correlators F ± k;s 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) 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 of (a/L), using the known value of the entries of the leading-order anomalous dimension matrix γ (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. (1);cont and X (1) g we have finally been able to compute the matrix γ (1);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:

NLO SF anomalous dimensions
• In our one-loop perturbative computation, we have obtained [X (1) O ] SF;lat for both c sw = 0 and c sw = 1 values. The results for [X (1) O ] cont;lat are known for generic values of c sw . We have thus been able to compute [X (1) O ] SF;cont for both c sw = 0 and c sw = 1 in such a way to check its independence from c sw . • For the "+" operators, we have checked the independence of γ (1);SF from the reference scheme used (either the RI-MOM or the MS). This is a strong check of the calculations from the literature of the NLO anomalous dimensions γ (1);cont and one-loop matching coefficients [X (1) O ] cont;lat in both the RI-MOM and MS scheme.
The resulting values of γ (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

Renormalization group running in perturbation theory
In this section we will discuss the perturbative computation of the RG running factorŨ (μ) 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 nonperturbative determination of the anomalous dimensions for QCD with N f = 2, the analysis below will be performed for that case; the qualitative conclusions are independent of the precise value of N f . The scale will be fixed using the value (20) MeV quoted in [38]. At leading order in perturbation theory the running factor is given by U LO in Eq. (2.33). Beyond LO, the running factor is given by Eq. (2.34), where W (μ) satisfies Eq. (2.36). In the computation of W , the β and the γ 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 γ /n β " where n γ is the order used for the γ function while n β is the order used for the β function. We will consider the following cases: (i) "1/1", i.e. the LO approximation in which W ≡ 1; (ii) "2/2", in which both γ and β are taken at NLO; (iii) "2/3", in which β is taken at NNLO and γ at NLO; (iv) "+3/3", in which β is taken at NNLO and for the NNLO coefficient γ 2 we use a guesstimate given by γ 2 γ −1 1 = γ 1 γ −1 0 ; (v) "−3/3", in which β is taken at the NNLO and for the NNLO coefficient γ 2 we use a guesstimate given by 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 β/γ 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 γ 2 = 0) in the 2/3 case, J 1 and J 2 with γ 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 γ and β 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 = γ 2 = 0) in the 2/2 case, J 3 (with b 3 = γ 3 = γ 2 = 0) in the 2/3 case, J 3 (with b 3 = γ 3 = 0 and γ 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.
In the Fierz − case we have thus proceeded by introducing the next order, namely: and γ 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 ±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 (μ) Eq. (2.41), we approximate better and better the numerical solution of Eq. (2.36), 12 which can thus be considered the best approximation of the running factor given a fixed-order computation of the β and γ 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 γ and β have simple expressions in terms of g(μ) rather than in terms of μ, Eq. (2.36) is most easily solved by rewriting it in terms of the derivative with respect to the coupling, viz.
where W (g(μ)) ≡ W (μ). While both terms on the righthand side diverge as g → 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 at some very perturbative coupling g i (but still a significantly larger value than required by Eq. (2.37)), where we include exactly the same coefficients J n , n = 1, . . . 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 γ /β used for the numerical integration. 13 Note that indeed the numerical value of g i needs not be extremely small for this to make physical sense, e.g. for N f = 2 (which will be of particular interest to us) and at the Planck scale one has g 2 MS (M P ) ≈ 0.221 ↔ α MS s (M P ) ≈ 0.0176 and g 2 SF (M P ) differs with respect to g 2 MS (M P ) only on the third decimal digit.
In Fig. 6b we compare the results for the numerical integration of W (μ) when matched at g 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 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 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 β and γ functions) obtained through numerical integration. They are presented in Table 1 at the scale μ = 3 GeV. In alternative we quote also the results for the 2/2 case perturbatively expanded at 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 μ = 3 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, O(g 2 ) case) from the other 3 (respectively, 4) numerical cases.
The results for the LO running factor U LO (μ) Eq. (2.33) and the numerically integratedŨ (μ) running factors beyond LO (ii)-(v) described above are illustrated in Figs. 7, 8, 9, 10, 11, and 12 together with the 2/2 O(g 2 ) perturbative expansion, for the four doublets of operators and three different schemes (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 γ 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 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 ±3/3 and have thus very large, asymmetric errors. • For both 4-5 submatrices (Fierz + and −) the ratio γ 1 γ −1 0 turns out to have large matrix elements. As a consequence, our plausibility argument for the guesstimates γ 2 γ −1 1 = ±γ 1 γ −1 0 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 MS and SF schemes, already for the 2/3 numerical solution (for the 2/2 O(g 2 ) perturbative expansion the situation is much worse). This obviously poses serious doubts on all the computations of F = 2 matrix elements beyond the Standard Model which uses perturbative running (in all cases through the 2/2 O(g 2 ) expansion) down to scales of 3 GeV or less.

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 f = 2 QCD will be presented in a separate publication [39]. 14 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 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 scalarscalar 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 (μ), 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 γ and β.
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 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 B = 2 amplitudes, should put a strong focus on non-perturbative renormalization. from chiral symmetry

In section 5.3 of [1] the authors derive an identity between the renormalization matrices for (Q
, 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 f |Q ± k |i , where Q ± k is a parity-even operator and |i, f are stable hadron states with the same, well-defined parity. Simple examples would be the matrix elements of F = 2 operators providing the hadronic contribution to K 0 -K 0 or B 0 -B 0 oscillation amplitudes (cf. Sect. 2). Bare matrix elements can be extracted from suitable three-point Euclidean correlation functions where O i, f are interpolating operators for the external states |i, f . If we perform a change of fermion variables of the form where ψ is a fermion field with N f flavor components and T is a traceless matrix acting on flavor space, this will induce a corresponding transformation Q ± k → Q k ± , O i, f → O i, f of the involved composite operators. If the regularized theory employed to define the path integral preserves exactly the SU(N f ) A axial chiral symmetry of the formal continuum theory, the equality 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 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 where ψ = (ψ 1 , ψ 2 , ψ 3 , ψ 4 ) T comprises the four, formally distinct flavors that enter Q ± k , Q ± k . Under this rotation, the ten operators of the basis in Eq. (2.1) transform as In the case of operators 1, 4, 5 the rotation is essentially trivial, in that it preserves Fierz (2 ↔ 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, 3 . In the latter expression we have written explicitly the renormalization scale μ. If we now use the RG evolution operators discussed in Sect. 2 to run Eq. (A.6) to another scale μ , one then has (recall that the continuum anomalous dimensions of Q + k and Q + k -respectively, Q − k and Q − k -are the same) which implies It is then immediate that the anomalous dimension matrices entering U ± are related as The correct interpretation of this identity is that, given an anomalous dimension matrix for, say, Q + 2,3 and Q + 2,3 , one can use Eq. (A.9) to construct a correct anomalous dimension matrix for Q − 2,3 and 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 γ − that follows from the condition imposed on Fierz − operators will lead to value of γ + 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 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 MS given in Appendix B, as well as (within uncertainties) the nonperturbative 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 γ SF 2 , 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 γ ± 0 , as can be readily checked in Eq. (2.24); therefore, the violation of the identity in e.g. SF schemes appears only at O(g 4 0 ) in perturbation theory.

Appendix B: NLO anomalous dimensions in continuum schemes
The two-loop anomalous dimension matrices in the RI-MOM scheme (in Landau gauge) [5,8] and MS scheme [8] are given by (the factor (4π) −4 has been omitted below to simplify the notation): Appendix C: Perturbative expansion of RG evolution for N 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 for N f = 30 and, more relevantly, for N f = 3 [46,47]. The reason is that, when Eq. (2.42) is written as a linear system, the 4 × 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 = 3 − N f ; if the NLO anomalous dimension matrix in the scheme under consideration is written as then one finds . (C.5) In the limit → 0 the element 23 of J ± 1 diverges; it is easy to check that the aforementioned 4 × 4 matrix, consistently, has determinant ∝ . A similar expansion of the matrices which is still divergent as → 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 (μ 2 , μ 1 ) = [Ũ (μ 2 )] −1Ũ (μ 1 ), on the other hand, is finite: the divergent part has the form M, and it is easy to check, using the explicit expression for U ± LO,0 (μ) and the identity U ± 15 The full expression for U (μ 2 , μ 1 ) in the → 0 limit still receives contributions from J 1,−1 , via the products with the O( ) terms in the expansion of U LO , which actually give rise to the only dependence of the expanded U (μ 2 , μ 1 ) on γ ± 1 . 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 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 f = 2 in detail. Indeed, applying exactly the same procedure for N f = 3 -i.e., solving Eq. (2.36) after having substituted the perturbative expressions for γ and β 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 f = 2 case we have employed W (μ 0 ) = 1 + g 2 (μ 0 )J 1 at some very high energy scale μ 0 . However, this can be replaced by the initial condition W (μ 0 ) = 1 at an even higher scale, thus again avoiding the appearance of any singularity; it turns out that the required value of g 2 (μ) 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 g 2 (μ) ∼ 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.
As a crosscheck of the robustness of our numerical approach we have computed explicitly the function W (μ) for N f = 3, using our numerical integration and W = 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 g 2 , to check whether the first coefficient J 1 is compatible within systematic fit errors (obtained by trying different polynomial orders up to O(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   (8) , 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 f . As a further crosscheck, we have also compared the result of computing the N f = 2 evolution with the two possible initial conditions. The outcome is that, if the value of the coupling at which W = 1 is sufficiently small, the two results are equal up to several significant figures down to values of the coupling g 2 2, where the hadronic regime is entered.
As showed in Sect. 3 the expected behavior of F( ) leads to the consideration of an asymptotic expansion of the form where the residue R n ( ) is expected to decrease faster as → ∞ than any of the terms in the sum. To determine the coefficients (α k , β k ) we minimize a quadratic form in the residues where F and ξ are the N − and (2n + 1)−column vectors (F( 1 ), . . . , F( )) T and (α 0 , α 1 , . . . , α n , β 1 , . . . , β n ) T , respectively, and f is the N × (2n + 1) matrix Again following [50], we have not introduced a matrix of weights in the definition of χ 2 . A necessary condition to minimize χ 2 is where we have assumed that the columns of f are linearly independent vectors (assuming 2n + 1 N ), and P is the  (6) 0.026133(6) −0.25536 (2) projector onto the subspace of R N generated by them. Equation (E.6) can be solved using the singular value decomposition of f , which has the form of where U is an N × (2n + 1) matrix such that With this decomposition one has Finally, the uncertainty in the result for ξ k can be modeled using error propagation as 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 χ 2 implies, in particular, that small values of might be given excessive weight. This problem has been dealt with by considering a range [ min , max ] with changing min . For this work the better convergence in results for (α k , β k ) was given by min = 16 and max = 46. The estimation of systematic uncertainty of the fitting procedure    (8) 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/ n+1 and log / n+1 with coefficients α n+1 and β n+1 , respectively. The systematic uncertainty of the finite part r 0 = α 0 is defined as the difference of the value of the parameter α 0 extracted by the two different fits. In the present work we have used n = 2 in the fit Ansatz for the O(a)-improved data, and n = 3 for unimproved ones (Tables 3, 4, 5, 6).