Non-Perturbative Renormalisation and Running of BSM Four-Quark Operators in $N_f = 2$ QCD

We perform a non-perturbative study of the scale-dependent renormalisation factors of a complete set of dimension-six four-fermion operators. The renormalisation-group (RG) running is determined in the continuum limit for a specific Schr\"dinger Functional (SF) renormalisation scheme in the framework of lattice QCD with two dynamical flavours ( $N_f = 2$ ). The theory is regularised on a lattice with a plaquette Wilson action and $\mathcal{O}(a)$-improved Wilson fermions. For one of these operators, the computation had been performed in ref. [1]; the present work completes the study for the rest of the operator basis, on the same simulations (configuration ensembles). The related weak matrix elements arise in several operator product expansions; in $\Delta F = 2$ transitions they contain the QCD long-distance effects, including contributions from beyond-Standard Model (BSM) processes. Some of these operators mix under renormalisation and their RG-running is governed by anomalous dimension matrices. In ref. [2] the RG formalism for the operator basis has been worked out in full generality and the anomalous dimension matrix has been calculated in NLO perturbation theory. Here the discussion is extended to the matrix step-scaling functions (matrix-SSFs), which are used in finite-size recursive techniques. We rely on these matrix-SSFs to obtain non-perturbative estimates of the operator anomalous dimensions for scales ranging from $\mathcal{O}(\Lambda_{\rm QCD})$ to $\mathcal{O}(M_W)$.


Introduction
In lattice QCD, the renormalisation of composite operators is an important step towards obtaining estimates of hadronic low-energy quantities in the continuum limit. Quark masses, decay constants, form factors, etc. are extracted from matrix elements of such operators; see ref. [3] for a recent review of lattice flavour phenomenology. Of interest to the present work is the class of dimension-six, four-fermion composite fields, arising in operator product expansions (OPE), in which the heavier quark degrees of freedom are integrated out. For ∆F = 2 and many ∆F = 1 transitions (F stands for flavour here), the resulting weak matrix elements of these operators govern long-distance QCD effects. They can be reliably evaluated by applying an intrinsically non-perturbative approach. Lattice QCD is our regularisation of choice which, by combining theoretical and computational methods, allows for an evaluation of these quantities with errors that can be reliably estimated and systematically improved.
Here we address the problem of calculating the renormalisation parameters and their renormalisation group (RG)-running for the operators defined in Eq. (2.1) below. We opt for the lattice regularisation consisting in the Wilson plaquette gauge action and the O(a)-improved Wilson quark action. We renormalise the bare operators in the Schrödinger Functional (SF) renormalisation scheme. This problem has first been studied with Wilson fermions for the relatively simple case of the multiplicatively renormalisable operators Q ± 1 of Eq. (2.1), both perturbatively [4] and non-perturbatively in the quenched approximation [5]. Subsequently results for Q ± 1 have also been obtained with N f = 2 dynamical sea quarks [1]. (An analogous study with quenched Neuberger fermions may be found in ref. [6].) Recently the perturbative calculations have been extended in ref. [2] for the rest of the operator basis Q ± k (k = 2, . . . , 5) of Eq. (2.1). The present is a companion paper of this work, complementing it by providing non-perturbative results for the operators Q ± k (k = 2, . . . , 5), computed in N f = 2 lattice QCD.
As stressed in refs. [2,7], these operators, treated here in full generality, become relevant for a number of interesting processes, once specific physical flavours are assigned to their fermion fields. For example, with ψ 1 = ψ 3 = s and ψ 2 = ψ 4 = d (cf. Eq. (2.2)), the weak matrix element K 0 |Q + 1 |K 0 comprises leading long-distance contributions in the effective Hamiltonian formalism for neutral K-meson oscillations in the Standard Model (SM). Allowing for beyond-Standard-Model (BSM) interactions introduces similar matrix elements of the remaining operators Q + 2 , . . . , Q + 5 . In some lattice regularisations, the corresponding bare matrix elements are expressed in terms of the operators Q + 1 , . . . , Q + 5 , with some important simplifications in their renormalisation properties [7][8][9]. In refs. [2] other flavour assignments are listed, leading to four-fermion operators related to the low-energy effects of ∆B = 2 transitions (B 0 -B 0 and B 0 s -B 0 s mixing) and to the ∆S = 1 effective weak Hamiltonian with an active charm quark. In the present work we will concentrate on the renormalisation and RG-runing of Q ± 2 , . . . , Q ± 5 . It is important to keep in mind that the sets {Q ± 2 , . . . , Q ± 5 } and {Q ± 2 , . . . , Q ± 5 } are parityeven and parity-odd components of operators with chiral structures (such as "left-left" or "left-right") which ensure their transformation under specific irreducible chiral representations. Chiral symmetry may be broken by the regularisation (e.g. lattice Wilson fermions) but it is recovered by the continuum theory. An important consequence is that our results, obtained for the continuum RG-evolution of the parity-odd bases {Q ± 2 , . . . , Q ± 5 }, are also valid for the parity-even ones {Q ± 2 , . . . , Q ± 5 }. The paper is organised as follows: In section 2 we list the operators we are studying and their basic renormalisation pattern. We also derive their RG-equations, and define the evolution matrices and the renormalisation-group invariant operators, which are scale-and scheme-independent quantities. This is an abbreviated version of section 2 of ref. [2]. The 1 interesting feature of the renormalisation pattern of operators Q ± k (k = 2, . . . , 5) is that they mix in pairs 1 . In the case of Q ± k (k = 2, . . . , 5), mixing is not an artefact of the lattice regularisation, as it also happens in schemes where all symmetries of the continuum target theory (QCD) are preserved; cf. ref. [7]. An important consequence of this property is that the RG-running of these operators is governed by anomalous-dimension and RG-evolution matrices, rather than scalar functions. The RG-evolution matrices are well known in NLO perturbation theory; cf. refs. [10,11]. Here, following ref. [2], we use them in closed form, suitable for non-perturbative evaluations. In section 3 we outline our strategy. First we define the SF renormalisation conditions for the operators Q ± k (k = 1, . . . , 5); again this is an abridged version of subsection 3.3 of ref. [2]. Next, we define in the SF scheme the matrix step-scaling functions (SSFs) as the RG-evolution matrices for a change of renormalisation scale by a fixed arbitrary factor; this factor is 2 in the present work. These are our basic lattice quantities, computed for a sequence of lattice spacings, at fixed renormalised gauge coupling. They have a well defined continuum limit, which is obtained by extrapolation, as explained in the same section. Repeating the calculation for a range of renormalisation scales (i.e. a range of renormalised couplings) and interpolating our data points, we finally have the SSFs as continuous polynomials of the gauge coupling, from which we obtain the anomalous dimension matrices, with NLO perturbation theory taking over only at O(M W ) scales. In section 4 we present our results. Besides the aforementioned SSFs, RG-evolution matrices and anomalous-dimension matrices, we also compute the renormalisation matrices for values of the gauge coupling corresponding to low-energy scales. These renormalisation factors are needed, in order to renormalise the corresponding bare matrix elements at these hadronic scales. The computation of the latter requires independent simulations on large physical lattices (of about 3-5 fm), which is beyond the scope of this work. Appendix A collects additional tests of the comparision between perturbative and nonperturbative RG evolution, including the specific renormalisation scale range, [2 GeV,3 GeV], considered in ref. [12]. Further details about one-loop cutoff effects in the SSFs are presented in Appendix B.

Renormalisation of four-quark operators
This section is an abridged version of sect. 2 of ref. [2], which we repeat here for completeness.

Renormalisation and mixing of four-quark operators
We recapitulate the main renormalisation properties of the four-fermion operators under study. These results have been obtained in full generality in ref. [7]. The absence of subtractions is elegantly implemented by using a formalism in which the operators consist of quark fields with four distinct flavours. A complete set of Lorentz-invariant operators is The operator subscripts obviously correspond to the labelling V → γ µ , A → γ µ γ 5 , S → 1, P → γ 5 , T → σ µν ,T → 1 2 ε µνρτ σ ρτ , with σ µν ≡ i 2 [γ µ , γ ν ]. Repeated Lorentz indices, such as γ µ γ µ and σ µν σ µν are summed over. In the above expression round parentheses indicate spin and colour traces and the subscripts 1, . . . , 4 of the fermion fields are flavour labels. Note that operators Q ± k are parity-even, and Q ± k are parity-odd. In the following we will assume a mass-independent renormalisation scheme. Renormalised operators can be written as (summations over l, m are implied), where the renormalisation matrices Z ± , Z ± are scaledependent and reabsorb logarithmic divergences, while ∆ ± , D ± are (possible) matrices of finite subtraction coefficients that only depend on the bare coupling. Throughout this work we use boldface symbols for the column vectors of four-fermion operator and the matrices which act on these vectors (e.g. Q, Q, Z, ∆, Z, D, etc.) while their elements are indicated with explicit indices (e.g. Q k , Q k , Z kl , ∆ kl , Z kl , D kl , etc.). We also introduce a simplification in our notation, by dropping the ± superscripts , wherever no ambiguity arises. This should not be a problem as the symmetric operator bases {Q + k } and{Q + k } (symmetric under flavour exchange 2 ↔ 4) never mix with the antisymmetric ones {Q − k } and {Q − k }, and thus equations are valid separately for each basis.
The renormalisation matrices have the generic structure Analogous expressions hold for Z and D. If chiral symmetry is preserved by the regularisation, both ∆ and D vanish. In the case of Wilson fermions, with chiral symmetry explicitly broken, we have ∆ = 0, whereas due to residual discrete flavour symmetries D = 0; this is the main result of ref. [7]. Therefore the left-left operators Q 1 = O VA+AV , which mediate Standard Model-allowed transitions, renormalise multiplicatively, while operators Q 2 , . . . , Q 5 , which appear as effective interactions in extensions of the Standard Model, mix in pairs: In conclusion, with Wilson fermions the parity-odd basis {Q k } renormalises in a pattern analogous to that of a chirally symmetric regularisation, while the parity-even one {Q k } has a more complicated renormalisation pattern due to the non-vanishing of ∆. We will henceforth concentrate on the non-perturbative renormalisation of the parity-odd basis {Q k } with Wilson fermions. 3

Renormalisation group equations
The scale dependence of renormalised quantities is governed by renormalisation group evolution. Denoting as µ the running momentum scale and µ the renormalisation scale where mass-independent renormalisation conditions are imposed, we have the following Callan-Symanzik equations for the gauge coupling and quark masses respectively: where f is a flavour label. The scheme mass-independence implies that the Callan-Symanzik function β and the mass anomalous dimension τ depend only on the coupling. Asymptotic perturbative expansions read Let us now turn to Euclidean correlation functions of gauge-invariant composite operators, of the form 2 with x = y j ∀j, y j = y k ∀j = k. For concreteness we have opted for correlation functions of the parity-odd operators Q k , which are the subject of the present work. Nevertheless, the results of this section apply to any set of operators that mix under renormalisation. The operators O l (l = 1, · · · , n) may be any convenient, multiplicatively renormalisable source field. For example they could be currents or densities (e.g. V µ (y), A µ (y), S(y) and/or P (y)), or Schrödinger functional sources at the time-boundaries. The latter will be explicitly discussed in Sect. 3. Renormalised correlation functions satisfy the system of Callan-Symanzik equations or, expanding the total derivative, where γ is a matrix of anomalous dimensions describing the mixing of {Q k } (cf. Eq. (2.12) below), andγ l is the anomalous dimension of O l (defined in a way analogous to Eq. (2.12)). A possible term arising from the running of the gauge parameter λ of the action is omitted here, for reasons explained in ref. [2]. A convenient shorthand notation for the anomalous dimension matrix of the operatorsQ k is thus The operator anomalous dimensions admit perturbative expansions of the form In standard fashion we can then derive This result implies that the block-diagonal form of the renormalisation matrices Z (and Z) of Eq. (2.4) induces the same block-diagonal structure for the anomalous dimension matrix γ. Thus the sums in Eqs. (2.12) and (2.14) simplify: for operatorQ 1 and its anomalous dimension γ 11 there is no summation; for operators Q 2 ,Q 3 summations run over indices 2 and 3 only, and similarly for the operator sub-basis Q 4 ,Q 5 .

Evolution matrices and renormalisation group invariants
In order to obtain a solution of Eq. (2.12) in closed form, it is convenient to introduce the renormalisation group evolution matrix U(µ 2 , µ 1 ) that evolves renormalised operators between scales 3 µ 1 and µ 2 < µ 1 : Substituting the above into Eq. (2.12) we obtain for the running of U(µ 2 , µ 1 ) with initial condition U(µ 1 , µ 1 ) = 1. Note that the r.h.s. is a matrix product. Following a standard procedure, the above expression can be converted into a Volterra-type integral equation and solved iteratively, viz.
where as usual the notation Texp denotes a Taylor expansion of the exponent, in which each term is an ordered (here g-ordered) product. Explicitly, for a generic matrix function M(x), we have

(2.18)
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 .
In the specific case of interest, M(g) = γ(g)/β(g), with γ(g) a matrix function and β(g) a real function. To leading order (LO) we have that M(g) = γ (0) /(b 0 g) and the independence of the matrix γ (0) from the coupling g simplifies eq. (2.18), so that the Texp becomes a standard exponential. One can then easily integrate the exponent in eq. (2.17) and obtain the LO approximation of the evolution matrix: When next-to-leading order corrections are included, the T-exponential becomes non-trivial. Further insight is gained upon realising that the associativity property of the evolution matrix U(µ 3 , µ 1 ) = U(µ 3 , µ 2 )U(µ 2 , µ 1 ) implies that it can actually be factorised in full generality as 20) and the matrixŨ(µ) can be expressed in terms of a matrix W(µ), defined through The matrix W can be interpreted as the piece of the evolution operator containing contributions beyond the leading perturbative order. Putting everything together, we see that and thus we make contact with the literature (see e.g. [10,11] W(µ) . (2.23) Expanding perturbatively we can check [2] that W is regular in the UV, and all the logarithmic divergences in the evolution operator are contained in U LO ; in particular, Rewriting Eq. (2.15) as and observing that the l.h.s. (respectively r.h.s.) is obviously independent of µ 1 (respectively µ 2 ), we conclude that these are scale-independent expressions. Thus we can define the vector of RGI operators aŝ where in the last step we use Eq. (2.24). 6

Schrödinger Functional renormalisation setup
In this section we introduce the finite volume Schrödinger Functional (SF) renormalisation schemes and the RG evolution matrix between scales separated by a fixed factor (i.e. the matrix-step-scaling function).

Renormalisation conditions
We first define Schrödinger Functional renormalisation schemes for the operator basis of Eq. (2.1). This section is an abridged version of sec. 3.3 of ref. [2]. We use the standard SF setup as described in [13], where the reader is referred for full details including unexplained notation.
We work with lattices of spatial extent L and time extent T ; here we opt for T = L. Source fields are made up of boundary quarks and antiquarks, where α, β are flavour indices, unprimed (primed) fields live at the x 0 = 0 (x 0 = T ) boundary, and Γ is a Dirac matrix. The boundary fields ζ,ζ are constrained to satisfy the conditions and similarly for primed fields. This implies that the Dirac matrices Γ must anticommute with γ 0 , otherwise the boundary operators O αβ [Γ] and O ′ αβ [Γ] vanish; thus Γ may be either γ 5 or γ k (k = 1, 2, 3).
Renormalisation conditions are imposed in the massless theory, in order to obtain a mass-independent scheme by construction. They are furthermore imposed on the parityodd four-quark operators {Q ± k } of Eq. (2.1), since working in the parity-even {Q ± k } sector would entail dealing with the extra mixing due to explicit chiral symmetry breaking with Wilson fermions, cf. Eq. (2.4). In order to obtain non-vanishing SF correlation functions, we 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 flavour indices, the minimal structure involves three boundary bilinear operators and the introduction of an extra, "spectator" flavour (labeled as number 5, keeping with the notation in Eq. (2.2)). We thus end up with correlation functions of the generic form where S s is one of the five source operators and similarly for S ′ s , which is defined with the boundary fields exchanged between time boundaries; e.g O 53 ↔ O ′ 53 etc. The constant η k is a sign that ensures F k;s (x 0 ) = G k;s (x 0 ) for all possible indices 4 ; it is easy to check that η 2 = −1, η s =2 = +1.We also use the two-point functions of boundary sources Finally, we define the ratios where α is an arbitrary real parameter. The geometry of F k;s , f 1 , and k 1 is illustrated in Fig. 1.
We can now impose Schrödinger functional renormalisation conditions on the ratio of correlation functions defined in Eq. (3.14), at fixed bare coupling g 0 , vanishing quark mass, and scale µ = 1/L. For the renormalisable multiplicative operators Q 1 we set For operators that mix in doublets, we impose 5 , (3.16) and similarly for Q 4,5 . The product of boundary-to-boundary correlators in the denominator of Eq. (3.14) cancels the renormalisation of the boundary operators in F k;s , and therefore Z jk;s 1 ,s 2 ,α only contains anomalous dimensions of four-fermion operators. Following [1,5,14], conditions are imposed on renormalisation functions evaluated at x 0 = T /2, and the phase that parameterises spatial boundary conditions on fermion fields is fixed to θ = 0.5. Together with the L = T geometry of our finite box, this fixes the renormalisation 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 restrict our choice to α = 0, 1, 3/2. One still has to check that the above renormalisation conditions are well-defined at treelevel. This is straightforward for Eq. (3.15), but not for Eq. (3.16): it is still possible that the matrix of ratios A has zero determinant at tree-level, rendering the system of equations for the renormalisation matrix ill-conditioned. This is indeed obviously the case for s 1 = s 2 , but the determinant vanishes also for other non-trivial choices of s 1 = s 2 . In practice, out of the ten possible schemes one is only left with six, viz. 6 (s 1 , (3.17) 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 ). Given the strong scheme dependence of the matrices γ (1);SF (cf. Eq. (2.13)), a criterion has been devised in ref. [2] in order to single out the scheme with the smallest NLO anomalous dimension. This consists in choosing the scheme with the smallest determinant and trace of the matrix 16π 2 γ (1);SF [γ (0) ] −1 for each non-trivial 2 × 2 anomalous dimension matrix. It turns out that the scheme defined by α = 3/2 and (s 1 , s 2 ) = (3, 5) satisfies these requirements in all cases (i.e. for the matrices related to (Q 2 , Q 3 ) and (Q 4 , Q 5 )). In the following we will present non-perturbative results for this scheme only 7 .

Matrix-step-scaling functions and non-perturbative computation of RGI operators
In order to trace the RG evolution non-perturbatively, we introduce matrix-step-scaling functions (matrix-SSFs), defined as 8 The above definition generalises the step-scaling functions (SSFs) defined for quark masses [14] and multiplicatively renormalisable four-fermion operators [5] such as Q ± 1 . Just like the anomalous dimension matrix γ, the matrix-SSF σ has a block-diagonal structure. So the above definition either refers to one of the two multiplicative operators Q ± 1 , or to one of the four pairs of operators that mix under renormalisation; i.e. (Q ± 2 , Q ± 3 ) or (Q ± 4 , Q ± 5 ). In the former cases σ is a real function, whereas in the latter cases it is a 2 × 2 matrix of real functions. Again in what follows the ± superscripts will be suppressed.
The advantage of working with step-scaling functions is that they can be computed on the lattice with all systematic uncertainties under control. More concretely, we define the lattice matrix-SSF Σ in a finite (L/a) 3 × (T /a) lattice; as repeatedly stated previously, in this work we set L = T . Working in the chiral limit, at a given bare coupling g 0 (i.e. at a given finite UV cutoff a −1 ) , Σ is defined as the following "ratio" of renormalisation matrices at two renormalisation scales µ = 1/L and µ/2 = 1/(2L): This quantity has a well defined continuum limit. For a sequence of lattice sizes L/a, we tune the bare coupling g 0 (a) (and thus the corresponding lattice spacing a) to a sequence of values which correspond to a constant renormalised squared couplingḡ 2 (1/L) = u. Keeping u fixed implies that the renormalisation scale µ = 1/L is also held fixed. It is then straightforward to check that Σ satisfies Thus, the computation of the renormalisation matrices Z at a fixed value of the renormalised squared coupling u and various values of the lattice sizes L/a and 2L/a, allows for a controlled extrapolation of the matrix-SSFs to the continuum limit. The strategy for obtaining non-perturbative estimates of RGI operators proceeds in standard fashion: We start from a low-energy scale µ had = 1/L max , implicitly defined bȳ g 2 (1/L max ) = u 0 . The SSF σ(u) for the coupling, defined as σ(ḡ 2 (1/L)) = g 2 (1/(2L)), is known for N f = 2 from ref. [15]. Thus we generate a sequence of squared couplings (u 1 , . . . , u N ) through the recursion σ −1 (u n−1 ) = u n , and compute recursively the matrix-SSFs (σ(u 1 ), . . . , σ(u N )) which correspond to a sequence of physical lattice lengths (inverse renormalisation scales) (L max /2, . . . , L max /2 N ). This is followed by the computation of is thought of as a high-energy scale, safely into the perturbative regime, and µ had ∼ O(Λ QCD ) as a low-energy scale, characteristic of hadronic physics. The RGI operators of Eq. (2.26) can finally be constructed as follows: In other words, once we know the column of renormalised operators Q(µ had ) at a hadronic scale from a standard computation on a lattice of "infinite" physical volume (which is beyond the scope of the present paper), we can combine it with the non-perturbative evolution matrix [U (µ had , µ pt ) (which is the result of this work) and the remaining µ pt -dependent factors at scale µ pt (known in NLO perturbation theory from ref. [2]), to obtain the RGI operators 9 .
All factors on the r.h.s. must be known in the same renormalisation scheme, which here is the SF. The scheme dependence should cancel in the product of the r.h.s., sinceQ is scheme-independent. In practice a residual dependence remains due to the fact that W(µ pt ) is only known in perturbation theory (typically to NLO). Finally we stress thatQ depends, through the operators Q(µ), on the values of the quark masses; of course the result also depends on the flavour content of the QCD model under scrutiny (i.e. N f ). We mentioned above that the matrix W(µ pt ) is known in NLO perturbation theory from ref. [2]. This statement requires a brief elucidation: W(µ pt ) is obtained by numerically integrating Eq. (2.23), using the NLO (2-loop) perturbative result for γ and the NNLO (3loop) perturbative result for β. In what follows this will be abbreviated as NLO-2/3PT. In line with ref. [2], also the present work devotes considerable effort to the investigation of the reliability of NLO-2/3PT at the scale µ pt .

Matrix-step-scaling functions and continuum extrapolations
We now turn to some practical considerations concerning the extrapolation of Σ(u, a/L) to the continuum limit a/L → 0, from which we obtain σ(u); cf. Eq. (3.20). We stress that although fermionic and gauge actions are Symanzik-improved by the presence of bulk and boundary counter-terms, correlation functions with dimension-six operators in the bulk of the lattice, such as those defined in Eqs. (3.4) and (3.14), are subject to linear discretisation errors. Their removal could be achieved in principle by the subtraction of dimension-7 counter-terms, but their coefficients are not easy to determine in practice. We therefore expect linear cutoff effects and consequently fit with the Ansatz In analogy to ref. [16], we explore the reliability of the above extrapolations with the help of the lowest-order perturbative expression for Σ ij , which includes O(ag 2 0 ) terms. In general the perturbative series for the operator renormalisation matrices has the form [16] where in the limit a/L → 0 the coefficients Z (l) are l-degree polynomials in ln(L/a) up to corrections of O(a/L). In particular the coefficient of the logarithmic divergence in Z (1) is given by the one-loop anomalous dimension γ (0) , and thus we parametrise Z (1) as with θ = 0.5 and T /L = 1. It is now easy to see that the one-loop perturbative expression for the matrix-SSF is given by In the continuum limit (a/L → 0 withḡ 2 = u fixed) we have The quantity contains all lattice artefacts at O(g 2 0 ). Results for δ k (L/a) are reported in Appendix B. The "subtracted" matrix-SSF, defined as also tends to σ in the continuum limit, but has the O(aḡ 2 ) effects removed. We will also use this quantity when studying the reliability of the linear continuum extrapolations below.

Perturbative expansion of matrix-step-scaling functions
Once the continuum matrix-SSF σ(u) has been computed for N discrete values of the renormalised couplingḡ 2 (1/L) = u, it is useful to interpolate the data so as to obtain σ(u) as a continuous function. This is done by fitting the N points by a suitably truncated polynomial With only a few (N ) points at our disposal, the fit stability is greatly facilitated by fixing the first two coefficients (matrices) r 1 and r 2 respectively to their LO and NLO perturbative values, leaving r 3 as the only free fit parameter. We will now derive the perturbative coefficients r 1 and r 2 .
Since the operator RG-running is coupled to that of the strong coupling, we also need the LO and NLO coefficients of its step-scaling function (SSF); i.e. (3.32) Given the strong coupling valueḡ 2 (1/L) = u at a renormalisation scale µ = 1/L, its SSF is defined as σ(u) =ḡ 2 (1/2L); cf. ref. [17]. Combining this definition with that of the Callan-Symanzik β-function of Eq. (2.5), we find that Plugging the NLO expansion of Eq. (2.7) in the above and taking Eq. (3.32) into account, we obtain the coefficients of the coupling SSF Matrix-SSFs for four-quark operators have been introduced in Eq. (3.18). In order to calculate the coefficients r 1 and r 2 of its perturbative expansion Eq. (3.31), we first write down the LO evolution matrix as Furthermore, the matrix W(µ) of Eq. (3.18) has the NLO perturbative expansion (cf. ref. [2] and references therein) from which the inverse matrix is readily obtained: We arrive at the last expression on the rhs by inserting the power-series expansion of σ(u) form Eq. (3.32). Substituting the various terms in Eq. (3.18) by the perturbative series (3.36), (3.37) and (3.38), we find From the first expression obtained for r 2 we see explicitly that O(u 2 ) corrections to W do not contribute (i.e. terms with J 2 are absent), in accordance with the fact that the O(u) term of W already contains all NLO contributions. The second expression for r 2 is obtained by using the property (cf. ref. [2] and references therein) Remarkably, the final result for r 2 is the exact analogue of the one found for operators that renormalise multiplicatively, cf. e.g. Eq. (6.6) in [4].

Non-perturbative computations
Our simulations are performed using the lattice regularisation of QCD consisting of the standard plaquette Wilson action for the gauge fields and the non-perturbatively O(a) improved Wilson action for N f = 2 dynamical fermions. The fermion action is Clover-improved with the Sheikoleslami-Wohlert (SW) coefficient c sw determined in [18]. The matrix-SSFs are computed at six different values of the SF renormalised coupling, corresponding to six physical lattice extensions L (i.e six values of the renormalisation scale µ). For each physical volume three different values of the lattice spacing a are simulated, corresponding to lattices with L/a = 6, 8, 12; this is achieved by tuning the bare coupling g 0 (a) so that the renormalised coupling (and thus L) is approximately fixed. At the same g 0 (a) we also generate configuration ensembles at twice the lattice volume; i.e. 2L/a = 12, 16, 24 respectively. We compute Z(g 0 , a/L) and Z(g 0 , a/(2L)) and thus Σ(g 2 0 , a/L); cf. Eq. (3.19). The gauge configuration ensembles used in the present work and the tuning of the lattice parameters (β, κ) are taken over from ref. [19] where all technical details concerning these dynamical fermion simulations are discussed. As pointed out in [19], the gauge configurations at the three weakest couplings have been produced using the one-loop perturbative estimate of c t [20], except for (L/a = 6, β = 7.5420) and (L/a = 8, β = 7.7206). For these two cases and for the three strongest couplings the two-loop value of c t [21] has been used.
Statistical errors are computed by blocking (binning) the measurements of each renormalisation parameter and calculating the bootstrap error on the binned averages. In order  to take their autocorrelation length into account, we determine the block-size for which the bootstrap error of a given renormalisation parameter reaches a plateau. This varies for each of the four matrix elements of a given 2 × 2 renormalisation matrix. We conservatively fix our preferred block-size to the maximum of all four cases, and estimate our statistical error accordingly. We crosscheck our results by also applying the Gamma method error analysis of ref. [22], and by varying the summation-window size. The results from the two methods agree within the (relevant) uncertainties. Numerical results for [Z(g 0 , a/L)] −1 and Z(g 0 , a/(2L)), computed from Eq. (3.16), are collected in Tabs. 7 and 8. The reason we prefer quoting the inverse of Z(g 0 , a/L) is that it is this quantity which is required for the computation of the matrix-SSFs; cf. Eq. (3.19).

Lattice computation of matrix-functions
We perform linear extrapolations in a/L of both Σ andΣ (cf. Eqs. (3.23) and (3.30)), so as to crosscheck the reliability of the continuum value σ(u). The extrapolation results can be found in Tabs. 1 and 2, as well as in Figs. 6,7,8, and 9. In most cases both extrapolations agree; at worst the agreement is within two standard deviations (e.g. in Fig. 6 the difference between off-diagonal elements of the matrices Σ andΣ is sizeable). We quote, as our best results, those obtained from linear extrapolations in a/L, involving all three data-points of the "subtracted" matrix-SSFs. We estimate the systematic error as the difference between the value of σ obtained by extrapolating Σ andΣ. This error is added in quadrature to the one from the fit.
Similar checks with another two definitions of "subtracted" matrix-SSFs, namely: which differ at O(u 2 ) have not revealed any substantial differences in the results.

RG running in the continuum
In order to compute the RG running of the operators in the continuum limit, matrix-SSFs have to be fit to the functional form shown in Eq. (3.31). Several fits have been tried out, with different orders in the polynomial expansion and r 2 either kept fixed to its perturbative value or allowed to be a free fit parameter. Fits with r 1 fixed by perturbation theory and r 2 the only free fit parameter do not describe the data well. This is understandable, as  deviations from LO are large for some matrix elements (for σ + 54 in particular) and knowledge of the NLO anomalous dimension γ (1) (and therefore r 2 ; cf. Eq. (3.40)) is necessary for a well-converging fit. It is however an encouraging crosscheck that the r 2 value returned by the fit is close to the perturbative prediction of Eq. (3.40). If, besides r 2 , we also include r 3 as a free fit parameter, the results have large errors. The best option turns out to be the one with the polynomial expansion of Eq. In the same Figures we also show the LO and NLO perturbative results, calculated from Eq. (3.31), truncated at O(u) and O(u 2 ) respectively. The comparison between the non-perturbative, the LO, and the NLO results provides a useful assessment of the reliability of the perturbative series. There is coincidence of all three curves at very small (perturbative) values of the squared gauge coupling u, but this is obviously guaranteed by the form of our fit function, as described above. At larger u-values one would ideally hope to see the NLO curves lying closer to the non-perturbative ones, compared to the LO curves. For σ + this is mostly the case, as shown in Fig. 2 55 curves are monotonically increasing, as opposed to the non-perturbative ones. In conclusion the overall picture in the renormalisation scheme under investigation is in accordance with our general expectations, although there are signs of slow or bad convergence of the perturbative results to the non-perturbative ones.
Once the matrix-SSFs are known as continuum functions of the renormalised coupling, we can obtain the RG-running matrix U(µ had , 2 n µ had ) = σ(u 1 ) . . . σ(u n ); cf. Eq. (3.21). We The matrixŨ(µ had ) does not depend on the higher-energy scale 2 n µ had , so the n-dependence on the rhs should in principle cancel out. We check this by computing the second line for varying n, using our non-perturbative result for U(µ had , 2 n µ had ) and the perturbative one forŨ(2 n µ had ). As explained in the comments following Eq. for The first error refers to the statistical uncertainty, while the second is the systematic one due to the use of NLO-2/3PT at the higher scale 2 n µ had . We estimate the systematic error as the difference between the final result, obtained with perturbation theory setting in at scale 2 8 µ had , and the one where perturbation theory sets in at 2 7 µ had (cf. Tabs. 3,4).
We note that systematic errors are almost negligible compared to statistical ones, the latter being the result of error propagation in the product of matrix-SSFs from µ had to 2 8 µ had . This however does not tell us much about the accuracy of NLO-2/3PT around the scale µ pt = 2 n µ had . We investigate this issue in Appendix A, where we compare σ(u n ), calculated in NLO-2/3PT and non-perturbatively. For several matrix elements of σ(u n ) we see that NLO-2/3PT is not precise enough, even at the largest scale we can reach (corresponding to n = 8).
We now play the inverse game, keeping fixed µ pt = 2 8 µ had and calculating It is computed for a fixed low-energy scale µ had and varying higher-scales 2 n µ had . For sufficiently large n, the results should not depend on the higher-energy scale. the computation thus described enforces the coincidence of our most perturbative point to the perturbative prediction, which we assume to describe accurately the running from µ pt ∼ O(M W ) to infinity. The discrepancies between perturbation theory and our results are evident at ever decreasing scales µ. These discrepancies are sometimes dramatic; e.g.
Finally, we compare the perturbative (NLO-2/3PT) to the non-perturbative RG evolution U(µ, µ * ) between scales µ and µ * , where µ * = 3.46 GeV is kept fixed and µ is varied in the range [0.43 GeV, 110 GeV]. The comparison is described Appendix A and confirms the unreliability of the perturbative computation of the RG running at scales of about 3 GeV.

Matching to hadronic observables with non-perturbatively O(a) improved Wilson fermions
Having computed the non-perturbative evolution matricesŨ(µ had ) as in Eq. (4.8), which provide the RG-running at the low energy scale µ had , we proceed to establish the connection between bare lattice operators and their RGI counterparts. Starting from the definition of Eq. (2.26), we write the RGI operator aŝ Z(g 2 0 , aµ had ) Q(g 2 0 ) .  Q is independent of any renormalisation scheme or scale; of course it is also independent of the regularisation. It is a product of several quantities: • The factors [g 2 (µ pt )/(4π)] − γ (0) 2b 0 and W(µ pt ) depend on a high-energy scale µ pt and are calculated in NLO perturbation theory. This was one of the main objectives of ref. [2].
• The running matrix U(µ pt , µ had ) is known between the high-energy scale µ pt and a low-energy scale µ had ; its non-perturbative computation for N f = 2 QCD is the main objective of the present work.
• The product of the last two factors Z(g 2 0 , aµ had ) Q(g 2 0 ) stands for the usual lattice computation of bare hadronic quantities and their renormalisation constants on large physical volumes and for several bare couplings, with the continuum limit taken though extrapolation.
Although the last item in the above list is beyond the scope of this paper, we have computed Z(g 2 0 , aµ had ) following [19], at three values of the lattice spacing, namely β = 6/g 2 0 = {5.20, 5.29, 5.40}, which are in the range commonly used for simulations of N f = 2 QCD in physically large volumes. The results are listed in Tabs. 5, 6. In order to interpolate to the target renormalized coupling u(µ had ) = 4.61, the data can be fitted with a polynomial. Our numerical studies reveal that additional values of β would be needed to improve the quality of the interpolation to the target value of the coupling.

Conclusions
In the present work we have studied the non-perturbative RG-running of the parity-odd, dimension-six, four-fermion operators Q ± 2 , . . . , Q ± 5 , defined in Eqs. (2.1) and (2.2). Assigning physical flavours to the generic fermion fields ψ 1 , . . . , ψ 4 , the above operators describe four-quark effective interactions for various physical processes at low energies. Under renormalisation, these operators mix in pairs, as discussed in Section 2. This mixing is not an artefact of the eventual loss of symmetry due to the (lattice) regularisation; rather it is a general property of operators belonging to the same representations of their symmetry groups. It follows that also the RG-running of each operator is governed by two anomalous dimensions, and the corresponding RG-equations are imposed on 2 × 2 evolution matrices. This makes the problem of RG-running more complicated than the cases of multiplicatively renormalised quantities, such as the quark masses or B K .
The novelty of the present work is that, using long-established finite-size scaling techniques and the Schrödinger Functional renormalisation conditions described in Section 3, we . The accuracy of our results for the diagonal matrix elements ranges from 3% to 5%. The accuracy on the determination of the non-diagonal matrix elements ranges from as high as 3% to as poor as 60%. Clearly there is room for improvement. In our next project concerning the renormalisation and RG-running of the same operators for QCD with three dynamical flavours, we plan to introduce several novelties, which ought to improve the precision of our results significantly. Perturbation theory is to be used for the RG-running for scales above µ pt ∼ O(M W ). In our SF scheme the perturbative results at our disposal are NNLO (3-loops) for the Callan-Symanzik β-function and NLO (2-loops) for the four-fermion operator anomalous dimensions. In Figs. 4 and 5 we see the presence of possibly relevant non-perturbative effects already at scales of about 3 GeV, where it is often assumed that beyond-LO perturbation theory converges well 10 . We have also performed some checks by computing the RG-evolution matrix from a generic scale to a scale of about 3 GeV and found some matrix elements where the NLO perturbative result significantly differs from the non-perturbative one (see Appendix A). This should serve as a warning for other non-perturbative approaches which assume that perturbation theory is convergent at such scales.
Finally, at a fixed hadronic scale and for three values of the bare gauge coupling, we have computed the renormalisation constants (again in 2×2 matrix form) of our four-fermion operators.
As a closing remark we wish to point out that the non-perturbative evolution matrices computed in this work describe not only the RG-running of the parity-odd operators Q ± 2 , . . . , Q ± 5 , but also that of their parity-even counterparts Q ± 2 , . . . , Q ± 5 . This is because evolution matrices are continuum quantities: in the continuum, each parity-odd operator combines with its parity-even counterpart to form an operator which transforms in a given chiral representation, both parts having consequently the same anomalous dimension matrices.
In the case, for instance, of ∆S = 2 transitions, we are dealing with operator matrix elements between two neutral K-meson states and therefore only the parity-even operators (Q + 1 in the SM and Q + 2,··· ,5 for BSM) contribute. Our results for the continuum RG-evolution, obtained for the parity-odd basis, can be used in this case. The renormalization of the bare operators, however, depends on the details of the lattice action. If the lattice regularisation respects chiral symmetry (e.g. lattice QCD with Ginsparg-Wilson fermions), then the parityeven and parity-odd parts of a given basis of chiral operators renormalise with the same renormalisation constants. Consequently they also have the same matrix-SSFs and evolution matrices. All results obtained for the parity-odd operators Q ± 2 , . . . , Q ± 5 are then also valid for the Q ± 2 , . . . , Q ± 5 , without further ado. Things are somewhat more complicated if the regularisation breaks chiral symmetry (e.g. lattice QCD with Wilson fermions). Then parity-even and parity-odd operators again have the same anomalous dimensions, as these are continuum quantities, but the "ratio" of their renormalisation matrices Z −1 Z is a finite (scale-independent) matrix which is a function of the bare gauge coupling; it becomes the unit matrix in the continuum limit. This "ratio" is fixed by lattice Ward identities, as discussed for example in ref. [7]. So the subtlety here is that once the renormalisation condition has been fixed for say, the parity-odd operator bases at a value g 2 0 of the squared gauge coupling, the condition for the parity-even coun-terparts is also fixed through Z −1 Z . Consequently, renormalisation matrix "ratios" like Z g 2 0 , a 2L Z g 2 0 , a L −1 are equal to their parity-even counterparts Z g 2 0 , a 2L Z g 2 0 , a L −1 .
Thus matrix-SSFs Σ(g 2 0 , a/L) and evolution matrices are the same for parity-odd and parityeven cases; cf. Eq. (3.19). But if we wish to use the evolution matrices of the present work also for the parity-even operators, we must ensure that these are renormalised in the "same" SF scheme employed for their parity-odd counterparts. This is ensured by writing the RGI parity-even operator column (in analogy to Eq. (4.9)) as: where Q sub ≡ (1 + ∆)Q is the "subtracted" bare operator, as suggested by Eq. (2.3). The term in square brackets of the last expression is the renormalised parity-even operator ZQ. It is computed however in a way that ensures that the bare operator Q(g 2 0 ) is renormalised in our SF scheme: the SF renormalisation parameter Z(g 2 0 , aµ had ) (which removes the logarithmic divergences) is multiplied by the scheme-independent, scale-independent "ratio" Z −1 Z .
Clearly, the procedure sketched above for the renormalisation of parity-even operators is fairly cumbersome. It is also prone to enhanced statistical uncertainties, as it involves subtracted operators Q sub with non-zero ∆. Fortunately, there is a way to circumvent the problem: it is well known that, using chiral (axial) transformations of the quark fields, we can obtain continuum correlation functions of specific parity-even composite operators in terms of bare correlation functions of parity-odd operators of the same chiral multiplet, regularised with twisted-mass (tmQCD) Wilson fermions [23]. The prototype example is the one expressing renormalized correlation functions of the axial current in terms of bare twistedmass Wilson-fermion correlation functions of the properly renormalised vector current. The situation is more complicated with four-fermion operators: in ref. [8] it was shown that such chiral rotations do indeed relate parity-even to parity-odd 4-fermion operators, but the resulting tmQCD Wilson-fermion determinant is not real, and thus unsuitable for numerical simulations. This problem is circumvented by working with a lattice theory with sea-and valence-quarks regularised with different lattice actions [9]. The valence action is the socalled Osterwalder-Seiler [24] variety of tmQCD, with valence twisted-mass fermion fields suitably chosen so as to enable the mapping of correlation functions involving parity-even operators {Q k } to those of the parity-odd basis {Q k }. The sea-quark action may be any tenable lattice fermion action. While the price to pay is the loss of unitarity at finite values of the lattice spacing, this is, however, outweighed by the advantage of vanishing finite subtractions (D = 0 in eq. (2.3)).

28
Appendix A Non-perturbative vs perturbative behaviour of the RG evolution In analogy to Appendix C of ref. [12] we construct the quantity: Once againŨ(2 n µ had ),Ũ(2 n+1 µ had ), and U(2 n µ had , 2 n+1 µ had ) are perturbative quantities known in NLO-2/3PT , while σ(u n+1 ) is a single non-perturbative matrix-SSF. In the last line of Eq. (A.1), the product [U(2 n µ had , 2 n+1 µ had ) −1 σ(u n +1 )] is the ratio of the nonperturbative over the perturbative RG evolution between scales 2 n µ had and 2 n+1 µ had . If perturbation theory were reliable at these high scales, D(n) would vanish at large n. The results for the D(n) matrix elements are shown in Figs. 10, 11. At the largest n values some of them are compatible with 0 while others are not. The latter case signals that due to large anomalous dimensions, NLO-2/3PT performs poorly even at scales as high as 2 n µ had and 2 n+1 µ had . Moreover, in Appendix C of ref. [12] the non-perturbative RG evolution U(µ, µ * ) between scales µ and µ * has been compared to the result from NLO-2/3PT. In ref. [12], µ is kept fixed to 2 GeV while µ * is varied in the range [2 GeV, 3 GeV]. We perform a similar study by fixing the reference scale µ * = 3.46 GeV = 2 3 µ had , corresponding to the squared coupling u 3 . This is the scale closest to the interval [2 GeV, 3 GeV] of ref. [12], for which we have directly computed the matrix-SSFs non-perturbatively. The scale µ is varied in the range [0.43 GeV, 110 GeV]. We compute U(µ, µ * ) in the following way: U(µ, µ * ) can be evaluated in a purely non-pertubative way for integer n 1 = log 2 (µ pt /µ) and n 2 = log 2 (µ pt /µ * ). The results are presented Figs. 12, 13. Very dominant non-perturbative effects are clearly visible for the elements (2, 2), (4,5), (5,4) and (5,5) of the operators Given the large deviation from the NLO-2/3PT running already seen in Fig. 5, these results are not surprising and simply confirm the non-reliability of the NLO-2/3PT computation of the RG running at scales around 3 GeV. Notice that the scale interval where this comparison has been performed in ref. [12] is completely contained in our plots between the third and the fourth point which correspond to scales of 1.73 GeV and 3.46 GeV. We remind the reader that a direct comparison between our results and those of ref. [12] is meaningless, the crucial differences being, among many others, the renormalisation scheme and the N f -value.    23 -0.14 -0.13 -0.12 -0.11 -0.10 -0.09 -0.08 -0.14 -0.

Appendix B One-loop cutoff effects in the step scaling function
In Tab. 9 we gather numerical values for δ k (L/a), defined in Eq. (3.29). We have calculated this quantity for a fermionic action with (c sw = 1) and without (c sw = 0) a Clover term. These results are also displayed in Figs. 14, 15, 16, and 17 (the target scheme α = 3/2, (s 1 , s 2 ) = (3, 5) is plotted with a blue triangle). Notice that the element (3, 2) of δ k is independent from α due to an accidental cancellation. This is why all data-points in the corresponding figures are not in colour. As expected, the Clover term has an important effect on the discretisation errors, which are significantly reduced when c sw = 1. The observed O(ag 2 0 ) discretisation effects in Figs. 14 and 15 are only due to the unimproved operators, the action being tree-level improved.