Nonperturbative Renormalization in Lattice QCD with three Flavors of Clover Fermions: Using Periodic and Open Boundary Conditions

We present the nonperturbative computation of renormalization factors in the RI'-(S)MOM schemes for the QCD gauge field ensembles generated by the CLS (coordinated lattice simulations) effort with three flavors of nonperturbatively improved Wilson (clover) quarks. We use ensembles with the standard (anti-)periodic boundary conditions in the time direction as well as gauge field configurations with open boundary conditions. Besides flavor-nonsinglet quark-antiquark operators with up to two derivatives we also consider three-quark operators with up to one derivative. For the RI'-SMOM scheme results we make use of the recently calculated three-loop conversion factors to the modified minimal subtraction scheme. The present version of the paper contains an Addendum with additional analytical expressions and updated results.


I. INTRODUCTION
Matrix elements of local operators between hadron states contain valuable information on hadron structure.For example, decay constants, moments of parton distributions and light-cone distribution amplitudes can be written in this form.In order to compute such quantities from QCD, nonperturbative methods are required.One possibility is lattice QCD: After a Wick rotation from Minkowski to Euclidean space the ensuing Euclidean functional integral of QCD is regularized with the help of a lattice discretization such that a numerical evaluation through Monte Carlo simulations becomes possible.
Hadronic matrix elements of local operators in general are ultraviolet (UV) divergent quantities.These are regularized on the lattice and their dependence on the cutoff is then treated by the renormalization procedure.In a perturbative continuum calculation, such as in dimensional regularization, the UV divergences are removed order by order in the expansion in the coupling constant by imposing a given condition on how the divergent part is subtracted (MS schemes).In lattice Monte Carlo simulations, expectation values and matrix elements of bare operators are determined nonperturbatively.Such quantities are however not interesting per se, but only as input for further computations in particle phenomenology.A calculation of the renormalization constants, e.g., in lattice perturbation theory, is therefore required to absorb the divergences of bare operators that appear as the lattice spacing a is sent to zero, and to convert the lattice results to the continuum schemes commonly adopted in phenomenology.
Unfortunately, the convergence of lattice perturbation theory is rather slow.Many strategies have been developed to address this limitation in the last decades.An alternative possibility is to employ an intermediate scheme which is applicable both in the continuum and on the lattice.The RI ′ -MOM scheme [1] is a Regularization Independent scheme where the renormalization condition is imposed on amputated vertex functions in momentum space at a renormalization scale determined by the external momenta.The renormalization constants in the RI ′ -MOM scheme can be calculated nonperturbatively on the lattice just by computing the expectation value of the vertex function on an ensemble of gauge-fixed configurations.The conversion (or matching) from the RI ′ -MOM to the MS schemes is then performed straightforwardly by a perturbative calculation in the continuum with dimensional regularization.
The continuum limit a → 0 is achieved by extrapolating results obtained at many nonvanishing lattice spacings a.As the lattice spacing is reduced, standard QCD simulation algorithms suffer from critical slowing down, meaning that an increasing Hybrid Monte Carlo (HMC) simulation time is required to fully sample the configuration space.If gauge (and fermion) fields fulfill (anti)periodic boundary conditions in all directions, the simulation will eventually get stuck in a fixed topological sector and ergodicity is lost.Open boundary conditions in the time direction have been proposed to solve this topological freezing problem at small lattice spacings [2].The modification of the space-time manifold allows the topological charge Q top to differ from an integer value, even in the continuum limit, thereby allowing for small continuous fluctuations of Q top during the HMC trajectory.This approach has been adopted within the CLS (coordinated lattice simulations) effort [3][4][5].
The breaking of translational invariance in the time direction is a drawback of open boundary conditions.Space-time correlators corresponding to physical particles can still be safely measured in the bulk of the lattice, far away from the boundaries.However, it is not a priori clear how the same procedure can be applied to vertex functions and correlators in momentum space.Here we show that this is indeed possible and present our approach to the computation of the renormalization constants in the RI ′ -(S)MOM schemes on n f = 2 + 1 CLS ensembles [3] with degenerate light and strange quark masses.
We consider a large variety of nonsinglet quark-antiquark operators with up to two derivatives as well as three-quark operators with up to one derivative employing lattices with lattice spacings down to about 0.04 fm.Treating the quarkantiquark operators in the RI ′ -SMOM scheme, we make use of the recently calculated three-loop conversion factors [6][7][8] in order to reduce the systematic error due to the truncation of the perturbative expansion.
With the help of Ward identities, the renormalization and improvement of the vector current is studied in Ref. [9] within the CLS setup.On a subset of the CLS ensembles with periodic boundary conditions some renormalization factors have been evaluated for a restricted range of lattice spacings in Ref. [10] using the RI ′ -MOM scheme.Utilizing Schrödinger functional techniques, the ALPHA collaboration has computed several (ratios of) renormalization factors within the CLS setup, see, e.g., Refs.[11][12][13][14][15][16].
The paper is organized as follows.In the next section we describe the quark-antiquark operators studied, the employed renormalization schemes as well as their numerical implementation on lattices with periodic and open boundary conditions.Section III discusses the renormalization of threequark operators.In Sec.IV we explain briefly the perturbative subtraction of lattice artifacts.Section V is devoted to the chiral extrapolation required to obtain a mass-independent renormalization procedure.Our conventions regarding continuum perturbation theory are collected in the following section.In Sec.VII we discuss the dependence of our results on the renormalization scale.The next section details the two basic methods that we apply to extract our final numbers.Some of these results are discussed in Sec.IX.Section X contains a short summary.Technical details and results from perturbation theory that are used in our computations are given in the Appendices A -G. Tables of our results can be found in Appendix H.The mixing matrices of the quark-antiquark and the three-quark operators are given in ancillary files.

II. QUARK-ANTIQUARK OPERATORS A. Multiplets of quark-antiquark operators
The main focus of our calculations is the renormalization of quark-antiquark operators involving up to two covariant derivatives acting on the quark fields.Three-quark operators, relevant for the computation of moments of baryon distribution amplitudes, can be treated analogously and will be discussed in Sec.III.
The elementary building blocks of local quark-antiquark operators are of the form where f , f ′ , . . .are flavor indices, α, β, . . .are spinor indices and Γ denotes a Dirac matrix.Color indices are suppressed and we define On the lattice the covariant derivatives are replaced by the standard (symmetric) discretized versions.We restrict ourselves to the flavornonsinglet case choosing, e.g., f ̸ = f ′ .Flavor-singlet operators could be constructed by summing over f = f ′ .
In order to retain as much of the continuum symmetry as possible on the lattice we consider multiplets of operators which transform irreducibly under the hypercubic group H(4) and enjoy a definite charge conjugation parity.Since the constraints imposed by space-time symmetry are less stringent on the lattice than in the continuum the possibilities for mixing increase.The choice of the operator multiplets is hence guided by the desire to avoid mixing as far as possible, especially mixing with operators of lower dimension.
When nonforward matrix elements of operators with derivatives are considered, either in the renormalization procedure or between hadronic states, we must be prepared to find mixing with so-called total-derivative operators, i.e., operators of the generic form This type of mixing occurs already in the continuum and is therefore unavoidable on the lattice.However, due to charge conjugation invariance, which is exact also on the lattice, we encounter this phenomenon only in the case of operators with two derivatives.The operator multiplets that we consider are compiled in Appendix A.

B. Renormalization schemes
In this section we describe our implementation of the RI ′ -MOM scheme [1] and the RI ′ -SMOM scheme [17].
Let O (i) m (x) (i = 1, 2, . . ., d, m = 1, 2, . . ., M ) denote M multiplets of local quark-antiquark operators which transform identically according to an irreducible, unitary, ddimensional representation of H (4). We call the unrenormalized, but (lattice-)regularized vertex functions (in Landau gauge) V (i) m (p, q), where p and q are the external quark momenta.The corresponding renormalized (in the MS scheme) vertex functions are denoted by V (i) m (p, q).The dependence of V (i) m on the renormalization scale µ is suppressed for brevity.Note that V (i) m as well as V (i) m carry spinor indices and are therefore to be considered as 4 × 4-matrices.(The color indices have been averaged over.)When all mixing multiplets are taken into account, we should have (up to power corrections in the lattice spacing) where Z is the matrix of renormalization and mixing coefficients and Z q is the wave function renormalization constant of the quark fields.
The renormalization factors depend on the renormalization scale µ as well as on the cutoff, i.e., the lattice spacing a.A third dimensionful quantity appearing in this connection is the asymptotic scale parameter Λ. Being dimensionless, the renormalization factors in a mass-independent renormalization scheme such as the RI ′ -(S)MOM scheme can only depend on two dimensionless combinations of these three quantities, e.g., on aµ and aΛ, or functions thereof, such as the bare lattice or the renormalized coupling constant.In the following, we shall suppress all arguments that are not needed in the respective context and write, e.g., Z(µ, a), Z(µ) or simply Z.
In the RI ′ -(S)MOM scheme we denote the matrix of renormalization and mixing coefficients by Ẑ and the quark field renormalization factor by Ẑq .The definition of Ẑq will be given below in Eq. (16).The renormalized vertex function is then written as With the lattice Born term B(i) m (p, q) corresponding to V (i) m (p, q), we obtain The RI ′ -(S)MOM scheme is now defined by requiring that for a given momentum geometry p = p, q = q the left-hand side of this equation coincides with the corresponding tree-level expression: In this way we arrive at our renormalization condition As the RI ′ -(S)MOM scheme should be mass independent, we have to impose this condition in the chiral limit.The sum over all members of the operator multiplets ensures that all lattice symmetries are preserved, if the individual operators are normalized such that they transform according to a unitary representation of H(4).
In the RI ′ -MOM scheme we choose while in the RI ′ -SMOM scheme we require p2 = q2 = (p − q) 2 = µ 2 , which may be achieved by taking Note that in our conventions the time component is the last component of the momenta.The choice of the momentum directions is to be considered as belonging to the definition of the renormalization scheme.It should be mentioned that the mixing with total-derivative operators can be taken into account only within the RI ′ -SMOM scheme, because these operators do not contribute in forward matrix elements.
If there is no mixing (M = 1), the formulas simplify.Omitting the superfluous multiplet indices m and m ′ in this case, we can write In the case of the (flavor-nonsinglet) local vector and axialvector currents one may want to modify the above renormalization conditions such that they are consistent with the respective Ward identities.Calling the vertex function of the local vector current V µ (with the flavor indices suppressed) V µ (p, q), we use the renormalization condition in the case of the RI ′ -MOM scheme.In the RI ′ -SMOM scheme we employ the renormalization condition [17] Ẑ−1 Similarly, we have for the axialvector current A µ with vertex function A µ (p, q) the renormalization condition in the RI ′ -MOM scheme and (15) in the RI ′ -SMOM scheme.
In the RI ′ -MOM scheme as well as in the RI ′ -SMOM scheme the wave function renormalization constant of the quark fields Ẑq is determined from the quark propagator S(p) according to with p2 = µ 2 .Other definitions of Ẑq have been proposed in the literature, mainly with the aim of reducing lattice artifacts, see, e.g., Refs.[17][18][19][20].
Using the lattice Born term instead of the continuum Born term in the renormalization condition (8) and proceeding analogously in the calculation of Ẑq ensures that Ẑ is the unit matrix in the free case.Note, however, that in many cases there is no difference between employing the lattice or the continuum Born terms for our choice of the momentum directions.
The renormalization matrix Ẑ leads from the bare lattice operators to renormalized operators in our regularization independent RI ′ -(S)MOM scheme.The matrix Z transforming the bare operators into renormalized operators in the MS scheme of dimensional regularization is then given by Z = C Ẑ, where the matrix C is calculated from Here m is the continuum Born term, and the factor C q is the analogue of C for the quark wave function renormalization constant, i.e., Z q = C q Ẑq .Note that C q and the conversion matrix C are completely determined from a calculation in continuum perturbation theory.If the renormalization of the vector or axialvector current is performed consistently with the Ward identities, cf.Eqs. ( 12) -( 15), the conversion factor C is equal to one in these cases.
One can avoid the use of the quark wave function renormalization constant by computing ratios Z/Z V with the help of the RI ′ -(S)MOM scheme and determining the renormalization constant Z V of the local vector current V µ by other methods, e.g., from matrix elements of V µ between hadronic states of given electric charge.
As already mentioned, all of the above calculations should be performed for massless quarks so that our renormalization schemes are mass independent.In practice, Monte Carlo simulations with our boundary conditions require nonvanishing quark masses.Therefore we compute Ẑ first at nonzero masses and perform an extrapolation to the chiral limit in the end, see Sec.V.In order to avoid unnecessary complications with this extrapolation we take the quark masses in the simulations to be flavor independent.

C. Numerical implementation
The CLS ensembles that we use are listed in Table I along with their most relevant properties.On the coarser lattices (β = 3.34, 3.40, 3.46, 3.55) we have ensembles with the standard boundary conditions at our disposal, i.e., periodic boundary conditions in all four directions for the gluons and periodic (antiperiodic) boundary conditions in space (time) for the quarks.In this case we use the term periodic boundary conditions as shorthand.The corresponding ensembles are labeled 'p' in Table I.On the finer lattices (β = 3.70, 3.85) only ensembles with open boundary conditions in time, labeled 'o', are available [3].It should be noted that the ensembles H101, U103, H200, and N202 are only used for the assessment of systematic uncertainties and not for the evaluation of our final results.
The lattice spacings have been determined from the Wilson flow time at the SU(3) symmetric point in lattice units t * 0 /a 2 , where the SU(3) symmetric point is defined by 12t * 0 m 2 π = 1.11.Equating t * 0 with the result µ * ref = (8t * 0 ) −1/2 ≈ 478 MeV of Ref. [21] we arrive at the values given in Table II.
Our calculations are based on correlation functions with external quark lines.Therefore we are forced to work in a fixed  gauge.As usual we choose the Landau gauge, because this gauge can be implemented on the lattice as well as in continuum perturbation theory.The required correlation functions are evaluated with the help of momentum sources.On the ensembles with periodic boundary conditions this procedure is straightforward to implement [22,23], and the relevant formulas are given in Sec.II C 1.In the case of open boundary conditions some modifications are necessary, which will be discussed in Sec.II C 2.

Periodic boundary conditions
For a given operator O(x) we start from the three-point function and the quark propagator where V denotes the (dimensionful) volume of the lattice.Since we use flavor-independent quark masses, flavor indices have been omitted.In the case of a flavor-nonsinglet operator O we have only quark-line connected contributions in the three-point function G, while for a flavor-singlet operator there would be an additional quark-line disconnected contribution.The quark propagator S always refers to a single flavor.
Note that due to translation invariance one of the sums in the above expression ( 18) is redundant.For example, one could restrict the sum over z to a single lattice point omitting at the same time a factor a 4 /V .However, the volume averaging connected with this sum suppresses statistical fluctuations very efficiently and is therefore highly advantageous.
Due to the invariance under global color transformations, which survives the Landau gauge fixing, both of these correlation functions are proportional to the unit matrix in color space.We assume that the color indices have been averaged over, as we already did in Sec.II B, so that G and S are 4 × 4 matrices.
The vertex function of the operator O is then constructed as The calculation of the correlation functions with the help of momentum sources proceeds as follows.If M (x, y) represents the fermion matrix on a given gauge field configuration, we compute the quark propagator Ŝ(x, y) with a momentum source by solving the lattice Dirac equation The quark propagator S(p) in momentum space is then evaluated as the gauge field average of Representing the quark-line connected part of G(p, q) is obtained as the gauge field average of The statistical errors are computed with the help of the (single elimination) jackknife procedure.As a relatively small number of configurations is sufficient for our purposes, we can choose them such that they are statistically independent to a good accuracy.Hence we refrain from binning.
If one imposes the standard boundary conditions on the quark fields, the possible spatial components of the quark momenta p on an N 3 s × N t lattice are integer multiples of 2π/(aN s ) while the time components are of the form (n + 1/2)2π/(aN t ) with n ∈ Z.Such momenta do not allow us to satisfy the condition (9) (or the condition (51) for three-quark operators below) exactly.Therefore we employ twisted boundary conditions [24] when solving the lattice Dirac equation (21).For twist τ in the spatial directions we get spatial momentum components p 1 = p 2 = p 3 = (n + τ /2)2π/(aN s ).In the RI ′ -MOM scheme the corresponding temporal twist is then chosen such that p 4 = p 1 = p 2 = p 3 .We employ the five values τ = 0.0, 0.4, 0.8, 1.2, 1.6.In the RI ′ -SMOM scheme with the momenta (10) twisted boundary conditions are not required.Nevertheless, we use the same twist values as in the RI ′ -MOM scheme, because a larger number of momenta appears to be beneficial in the analysis.
Besides the renormalization factors themselves one might also be interested in ratios of these factors, either in order to avoid the appearance of the quark wave function renormalization constant or because the physical quantity to be studied is a ratio of matrix elements of different operators and is hence renormalized by the ratio of the corresponding renormalization factors.The latter case appears for example when moments of hadron distribution amplitudes are computed [25][26][27][28][29].Such ratios of renormalization coefficients can be evaluated "directly" by computing them on the single jackknife ensembles and then analyzing the results as for individual renormalization factors.D4ψ on ensemble J500 summed over space as a function of time (in lattice units).The magnitude of the momentum is given by µ ≈ 2 GeV and a −1 ≈ 5 GeV.The blue and dotted grey lines are as in Fig. 1.

Open boundary conditions
On the ensembles with open boundary conditions we have to modify our computational strategy.In order to avoid instabilities in the inversion of the fermion matrix, i.e., in the evaluation of the quark propagators, the support of the momentum sources is no longer taken to be the whole N 3 s ×N t lattice, but replaced with a subvolume of size N 3 s × (N t − 2∆) situated symmetrically about the center of the lattice, thus keeping a distance of ∆ • a from the temporal boundaries.Moreover, we have to take into account that the invariance under translations in the time direction is broken.Still, in the bulk of the lattice, i.e., at sufficiently large distances from the temporal boundaries, the effects of this symmetry loss should be negligible.Therefore we restrict the sum over the operator position z in the three-point function (18) to an even smaller volume of size N 3 s × (N t − 2 ∆), again centered at the midpoint of the lattice.
In our calculations we choose ∆ = 4 and ∆ = 32.To motivate these choices we plot spatial sums of correlation functions against time on the ensemble J500 in Figs. 1 and 2. In Fig. 1 we do this for the quark propagator (19) and in Fig. 2 for the three-point function (18) of the operator ψγ 1 ↔ D 4 ψ.In both cases we show the real parts of all 16 components for the momentum geometry (9) of the RI ′ -MOM scheme with µ ≈ 2 GeV.The dotted vertical lines indicate the source volume, while the blue vertical lines limit the subvolume that will be used for the final summation.It should be noted that the "flatness" of the data in the bulk improves with increasing µ and that µ = 2 GeV is the lowest scale entering our final analysis.We remark that J500 corresponds to our smallest lattice spacing.Therefore, at the other β values, ∆ • a = 32a is bigger if expressed in physical units.
Also the implementation of the Landau gauge requires some modifications on lattices with open boundary conditions.On a given lattice gauge field configuration the Landau gauge is imposed by maximizing the functional with respect to the gauge transformation field Ω(x) acting on the link variables U µ (x) as The corresponding Lie-algebra valued field A µ (x) is given by the traceless antihermitian part of the gauge-fixed link field U µ (x), i.e., The maximization can be performed, for instance, by a local overrelaxation algorithm.We continue the minimization iterations until the deviation from the Landau gauge condition, is smaller than 10 −10 .In the case of open boundary conditions we somewhat modify this procedure: In the sum in Eq. ( 25) an additional factor 1/2 is attached to the spatial links living on the boundaries at x 4 = 0 and x 4 = aN t , in analogy to the modification of the standard plaquette gauge action [2].

Finite size effects
The finite volume of our lattices may distort our results.On the one hand, we should therefore consider lattices of different size with otherwise identical parameters and perform an infinite volume limit in the end.This would amount to a rather demanding procedure.On the other hand, renormalization is a short-distance phenomenon.Hence one may expect that the evaluation of renormalization factors is not severely affected by finite size effects, at least for not too small renormalization scales.
In the end we are restricted to the ensembles in Table I, which do not allow us to perform a systematic infinite volume limit.However, for two simulation points (β = 3.40, κ = 0.13675962 and β = 3.55, κ = 0.137) we do have ensembles with different spatial volumes at our disposal so that we can get at least some hints at the size of finite volume effects.This is illustrated in Fig. 3 for the v 2b operators given in Eq. (A7).We show Z RGI , as defined in Eq. (67) below, obtained on a spatial volume N 3 s = 32 3 at β = 3.55, κ = 0.137 with periodic boundary conditions (blue squares) as well as on two lattices with open boundary conditions, 32 3 (black circles) and 48 3 (red triangles).The amount of agreement between the blue squares and the black circles indicates the consistency between the results from open and periodic boundary conditions, while the comparison of the black circles with the red triangles gives an impression of the finite size effects.Note that only results with µ 2 ≥ 4 GeV 2 will enter our final analysis.

III. THREE-QUARK OPERATORS
In this section we describe briefly the renormalization of three-quark operators, which appear in the description of baryon distribution amplitudes [25,26].
Every local three-quark operator can be represented as a linear combination of the operators Here we use a multi-index notation for the covariant derivatives, Dl ≡ D λ1 • • • D λ l .Aiming again at a mass-independent renormalization scheme we assign the same mass to all flavors and eventually consider the chiral limit.In an abbreviated notation we write the above operators as Then we have for all permutations π in the symmetric group S 3 of three elements, where From these "elementary" operators we construct the operators of interest with the help of flavor structures F and spinor structures S according to where a sum over all (multi-)indices that appear twice is implied.
In the case of the quark-antiquark operators one has to distinguish between flavor-singlet and flavor-nonsinglet operators corresponding to the decomposition 3 ⊗ 3 = 1 ⊕ 8 under flavor SU(3).For our three-quark operators we have the decomposition 3 ⊗ 3 ⊗ 3 = 1 ⊕ 8 ⊕ 8 ⊕ 10.The flavorsinglet (flavor-decuplet) representation corresponds to the totally antisymmetric (totally symmetric or trivial) representation of S 3 .The two flavor octets, called mixed symmetric (MS) and mixed antisymmetric (MA), form a basis for the two-dimensional representation of S 3 .More explicitly, we have the singlet flavor structure F B,f1f2f3 s with decuplet flavor structures F B,f1f2f3 d with and the octet flavor structures F B,f1f2f3 ot , where the second subscript t takes the value t = 1 for MS and t = 2 for MA.
The spinor structures should be chosen to yield a flavorspinor structure that is totally symmetric under simultaneous permutations of the flavor, spinor and derivative indices f a , α a , la (a = 1, 2, 3).Furthermore, the operator multiplets should transform irreducibly under the spinorial hypercubic group H(4), which replaces the hypercubic group H(4) in the case of fermionic operators [30].The group H(4) has five irreducible spinorial representations: τ (The superscripts indicate the dimension of these representations.)Multiplets transforming according to these representations have been given in Ref. [31].Starting from these operators we construct multiplets of spinor structures which transform under S 3 identically to their flavor counterparts: Here m labels the different H(4) multiplets and i labels the different members of the multiplets.Then is indeed totally symmetric under simultaneous permutations of the flavor, spinor and derivative indices.An analogous statement holds for the singlet and the decuplet.The corresponding operators S l α Ψ f α ( l; x) with zero and one derivative have been given in Ref. [25] for generic flavors.For the reader's convenience they are collected in Appendix B. They are chosen such that they do not mix with operators of lower dimension.Moreover, there is no mixing between operators transforming according to nonequivalent representations of S 3 .
In the case of the flavor-octet operators we find for t ∈ {1, 2}.Therefore we can always work with the MA flavor structure (t = 2) and assume that the flavor-spinor structure factorizes into a flavor structure and a spinor structure as in Eq. (33).For the singlet and decuplet operators this factorization is trivially satisfied.Our renormalization procedure for the three-quark operators is similar to the RI ′ -SMOM scheme used in the case of the quark-antiquark operators.In particular, we compute the quark field renormalization factor Z q from the quark propagator as above, see Eq. (16).
We write the mixing operator multiplets for a fixed flavor structure F B,f o2 in the form and analogously with o2 replaced by o1, s or d.The corresponding vertex functions are given by The renormalized vertex functions take the form where The renormalization and mixing coefficients Ẑmm ′ are fixed by the renormalization condition which is analogous to Eq. ( 8) for quark-antiquark operators.
The superscript "Born" indicates the corresponding tree level expression (Born term).This is again taken with all lattice artifacts included.More explicitly our renormalization condition can be written as with and For singlet and decuplet one gets analogous equations where, of course, no sums over t appear.So we have The corresponding (matrices of) renormalization factors leading to operators renormalized in the MS scheme are constructed with the help of (matrices of) conversion factors calculated in continuum perturbation theory, where we use the particular version of the MS scheme introduced in Ref. [32].Due to the complexity of higher-loop calculations we had to limit ourselves to one-loop accuracy [33].Also the anomalous dimensions of our operators are in general only known to one loop, with the exception of the operators without derivatives, for which the anomalous dimensions have been calculated to three loops [34].
The evaluation of the correlation functions on lattices with periodic and open boundary conditions as well as the chiral extrapolation proceed in complete analogy to the case of quark-antiquark operators.For the external momenta we have taken employing twisted boundary conditions.

IV. PERTURBATIVE SUBTRACTION OF LATTICE ARTIFACTS
For larger values of the renormalization scale µ lattice artifacts will show up.Given the fact that for most operators Z has to diverge as a → 0, it is not immediately obvious what one should call lattice artifacts in the present context.In order to clarify this point it is useful to have a look at the calculation of Z in lattice perturbation theory.Evaluating the required correlation functions to one-loop order at vanishing quark mass yields results of the form where g is the bare coupling constant and F (0) = 0. What is usually quoted as the one-loop result from lattice perturbation theory is the above expression with all contributions that go to zero for a → 0 omitted, i.e., From this point of view the quantity F (a 2 µ 2 ) would be considered as a lattice artifact, which vanishes for fixed µ like a 2 (up to logarithms) as a → 0 [35].However, keeping such contributions (or part of them) should not change the continuum limit of the renormalized quantities.Nevertheless, it seems to be generally expected that suppressing the above sort of lattice artifacts in the renormalization factors would also reduce the lattice artifacts in the renormalized quantities.
One can realize such a suppression by calculating expressions for the lattice artifacts in (one-loop) lattice perturbation theory and subtracting these from the data [23,36].Using the above notation this means subtracting We have computed F (a 2 µ 2 ) for the quark-antiquark operators with less than two derivatives in the RI ′ -MOM and the RI ′ -SMOM schemes (see Appendix C for details).
Of course, the bare coupling g 2 in this expression could be replaced by other sensible definitions of the coupling constant such as the boosted coupling where 1  3 tr U denotes the average plaquette.In most cases, both choices of the coupling constant seem to reduce the discretization effects by about the same amount, though with opposite signs of the remaining lattice artifacts.
In the following we shall restrict ourselves to the straightforward case of the bare coupling.A scale-dependent choice of the coupling has been proposed in Ref. [10].

V. CHIRAL EXTRAPOLATION
Aiming at a mass-independent renormalization scheme we should use massless quarks.However, the boundary condi-tions that we employ in our Monte Carlo simulations require massive quarks.Hence we must extrapolate the results obtained in our simulations with three mass-degenerate quarks to the chiral limit, where all quark masses vanish.
At β = 3.34 there are ensembles with periodic boundary conditions for three different quark masses available, while for each of the next three β values (β = 3.40, 3.46, 3.55) we have ensembles with periodic boundary conditions for four different masses, see Table I.On the two finest lattices (β = 3.70, 3.85) we have only two different masses each at our disposal.All of these ensembles will be used for the chiral extrapolation.
One-loop lattice perturbation theory [18] suggests that the leading mass dependence is linear in the quark mass.Therefore, in a first approach, we extrapolate linearly in m 2 π .On the coarser lattices we observe a rather mild mass dependence of the vertex functions, especially for larger scales, so that we are confident that the linear extrapolations yield reliable results also on the two finest lattices.The chiral extrapolations are performed at fixed external momenta.As these depend on the lattice size (see Sec. II C 1) it is in some cases necessary to interpolate between the simulated momenta.This is done linearly in µ 2 .As an example we show the chiral extrapolation of the tensor density at β = 3.55 for three scales in Fig. 5.
On the coarser lattices we could perform additional quadratic (in m 2 π ) extrapolations in order to estimate the uncertainties connected with taking the chiral limit.However, this is not possible on the two finest lattices.In order to get nevertheless an impression of the importance of m 4 π terms in the full range of β we have used the following procedure as a second approach.
Since the renormalization factors computed in the massless theory suffice to renormalize also the vertex functions evaluated with nonvanishing quark masses, it follows that has a finite continuum limit a → 0.Here we have indicated the dependence of the renormalization matrix Ẑ on the renormalization scale µ, the lattice spacing a and the pion mass m π explicitly.Note that µ has to be kept fixed in physical units as a → 0. The existence of the continuum limit of Ẑ(µ, a, m π ) Ẑ(µ, a, 0) −1 motivates the following ansatz for the mass dependence: Note that the coefficients b 0 , b 1 , c 0 , and c 1 are of mass dimension −2, −1, −4, and −3, respectively.One-loop lattice perturbation theory [18] reveals contributions proportional to am q .Therefore we include terms linear in a to describe the lattice spacing dependence of Ẑ(µ, a, m π ) Ẑ(µ, a, 0) −1 .Since we have to work at the same value of µ for all lattice spacings, it is necessary to interpolate the scale dependence of the Ẑ data.For this purpose we use cubic splines in ln(a 2 µ 2 ).
If we consider M mixing operator multiplets, the quantities b parameters, if we use data for n a values of a.In the case of the quark-antiquark operators, where n a = 6, these 10M 2 parameters are fitted to 19M 2 data points corresponding to the 19 available combinations of a and m π .In the case of the three-quark operators we do not have data for β = 3.34 at our disposal, hence n a = 5, and 9M 2 parameters must be fitted to 16M 2 data points.For M > 2 the number of free parameters becomes prohibitively large, and we have to restrict ourselves to separate fits for each β.
In the following we shall call the first approach "local chiral extrapolation" because the extrapolation is performed independently for each β.The second approach will be referred to as a "global chiral extrapolation".The latter method addresses the limitation that we do not have the same coverage of pion masses at all the lattice spacings, while the former approach has been used in our previous publications.One finds rather good agreement between the two procedures, which improves for larger scales.An example is shown in Fig. 6.Nevertheless we consider the global chiral extrapolation to be more reliable, because in this case we expect the extrapolation on the two finest lattices to benefit from the information on the behavior at small masses, which is available only on the coarser lattices.

VI. INPUT FROM CONTINUUM PERTURBATION THEORY
Continuum perturbation theory is used for the calculation of the conversion factors.Moreover, in order to control the scale dependence of the renormalization factors we need perturbative expressions for W (µ, µ 0 ) = Z(µ)Z −1 (µ 0 ).In both cases the running coupling ḡ(µ) is required.
The scale dependence of ḡ(µ) is controlled by the β function where the derivative is to be taken at fixed bare parameters and fixed cutoff.The perturbative expansion of β(ḡ) can be written as The values of the coefficients β 0 -β 4 in the MS scheme are listed in Eqs.(D1) -(D5).Integrating Eq. ( 58) one obtains with the Λ parameter appearing as an integration constant.The running of the Z matrices is governed by the anomalous dimension matrix whose perturbative expansion reads (62) Again, the bare parameters and the cutoff are kept constant when taking the derivative.
We define the quark field renormalization constant Z q so that the renormalized quark propagator is obtained from the bare propagator by multiplication with Z q , i.e., we use the continuum analogue of Eq. ( 16).For the anomalous dimension of the quark field we adopt the definition Collections of coefficients γ i in the MS scheme can be found in Appendix D for our quark-antiquark operators as well as for the quark field and in Appendix E for our three-quark operators.
Taking the derivative with respect to the running coupling ḡ(µ), we get This system of differential equations can formally be solved in the form In the general case of M mixing multiplets of operators this expression may be difficult to evaluate, because the M × M matrices γ(g) do not necessarily commute for different values of the coupling g.If there is no mixing (M = 1) the formula simplifies to In analogy to the definition of the Λ parameter we can define the so-called RGI (renormalization group invariant) Z, which is independent of scale and scheme: Of course, the independence of scale and scheme is strictly realized only if the β function and the anomalous dimension γ are exactly known (including nonperturbative contributions).In practical applications this is usually not the case, except for the vector current and the nonsinglet axialvector current, whose anomalous dimensions are known to vanish exactly.In all the other cases one must be prepared to encounter violations of the scale and scheme independence due to truncation errors in the perturbative approximations.
With the help of the methods described in Secs.II B and II C we can compute Ẑ, the renormalization matrix leading from the bare lattice operators to operators renormalized in the RI ′ -(S)MOM scheme, for some range of renormalization scales µ.In order to evaluate the matrix Z transforming the bare operators into MS operators we need in addition the conversion matrix (or conversion factor if there is no mixing) C, cf.Eq. ( 17) in the case of quark-antiquark operators.This is calculated in continuum perturbation theory leading to results of the form For selected quark-antiquark operators, explicit expressions for the coefficients c 1 , c 2 , . . .are compiled in Appendix F. For three-quark operators, see Appendix G.
With the β function and the anomalous dimension function γ given to some order in perturbation theory, we calculate the integrals in ( 60) and ( 66) analytically if possible or else by numerical integration.In the case of two mixing multiplets with a nondiagonal matrix of anomalous dimensions, the integral representation (65) of the scale dependence of Z is not very helpful.Instead of trying to evaluate the expression (65) we solve the system of differential equations (64) with the help of a power series expansion in ḡ following Ref.[37].
In principle, one has a lot of freedom in choosing the renormalization scheme for the anomalous dimension, the coupling constant in the conversion factor etc., cf.Ref. [23].If there were no truncation errors, Z RGI would be independent of all these choices.In practice, we consider only the MS scheme.In order to obtain an estimate of the truncation errors we vary the number of loops taken into account.

VII. LOOKING FOR A WINDOW
At this point, the results for the Z factors (extrapolated to the chiral limit and converted to the MS scheme) suffer from two problems: the truncation errors of the perturbative expressions for the γ and β functions and the conversion matrix (factor) C on the one side and the lattice artifacts on the other side.The former grow as µ becomes smaller, while the latter increase as µ gets larger.Ideally, one would like to find a window, i.e., a µ interval where both errors are negligible.This would require for a lattice of linear extent L. Then lattice artifacts would be negligible and the scale dependence could be described by low-order continuum perturbation theory.In such a window the results for Z RGI would be independent of the scale µ, i.e., when plotted against µ they would show a plateau, and the final value for Z RGI could be read off at a value µ within this window.
Unfortunately, such an ideal situation is hard to achieve at the present lattice spacings.Although in recent years there has been considerable progress in the perturbative evaluation of anomalous dimensions and conversion factors, the convergence of the perturbative expansions, in particular for the conversion factors, is not always satisfactory.Therefore the results for Z RGI will usually not be independent of µ even for values µ ≈ 2 GeV.
For increasing values of µ lattice artifacts will become larger.Nevertheless, the continuum limit of renormalized quantities remains in principle well-defined also for higher renormalization scales.Therefore one can simply use Z evaluated at some convenient, fixed scale µ and hope that the continuum limit is under control.However, one may expect that suppressing the lattice artifacts in the Z factors would also make the behavior of the renormalized quantities for a → 0 more benign.
One possibility to realize such a suppression consists in fitting the Z RGI data with a suitable ansatz for the lattice artifacts.Since the values of p 2 = a 2 µ 2 actually appearing in the data cover only a finite interval which does not extend down to 0, a polynomial in p 2 would be a reasonable choice.
Alternatively one can calculate expressions for the lattice artifacts in lattice perturbation theory and subtract these from the data as explained in Sec.IV.While this procedure reduces the scale dependence of the Z RGI data already substantially, some lattice artifacts will in general remain, which can still be fitted.
Lattice artifacts can be studied in a much clearer way for quantities that have a decent continuum limit.This is not the case for the Z factors themselves, but for ratios Z(µ)Z −1 (µ 0 ) of renormalization factors (matrices) at different renormalization scales [38].Therefore these ratios offer a possibility to investigate the impact of discretization effects in our calculations.After performing a continuum limit we can see how well the scale dependence is described by continuum perturbation theory.Such investigations have already been performed for various operators in a number of different settings, see, e.g., Refs.[39][40][41].In this way we should be able, at least for these quantities, to disentangle lattice artifacts and truncation errors, which otherwise are hard to separate unambiguously.
We show Z(µ)Z −1 (µ 0 ) for the tensor density in Fig. 7 comparing results obtained with conversion factors at different orders in perturbation theory.The statistical errors of the ratios have been calculated from the errors of numerator and denominator by means of error propagation.For the continuum extrapolation we have employed a second-order polynomial in a 2 .The agreement with the perturbative scale dependence shown by the curves improves significantly as we go from one-loop to three-loop conversion factors.In general the perturbative series for the conversion factors is less well behaved than the expansions of the anomalous dimensions, which tend to converge quite fast.When plotting Z RGI against µ we must therefore be prepared to find visible deviations from a plateau also for moderate values of µ, where lattice artifacts should be small.
If the one-loop lattice artifacts are not subtracted, the agree- The nonperturbative results at finite β have been computed in the RI ′ -SMOM scheme, lattice artifacts have been subtracted perturbatively, and a global chiral extrapolation has been performed.Finally, the values have been converted to the MS scheme using one-loop, two-loop and three-loop (from top to bottom) perturbation theory and extrapolated to the continuum limit (β = ∞).The curves show the behavior calculated in three-loop perturbation theory in the MS scheme.ment with the perturbative running is less satisfactory, even if the three-loop conversion factor is used.For the tensor density this is demonstrated in Fig. 8, which should be compared with the lowest panel in Fig. 7.
Particularly good agreement with the perturbative running is observed for the pseudoscalar density as Fig. 9 shows.In this case we can even use the five-loop expression for the anomalous dimension (see Appendix D), but the difference with the three-loop results would be hardly visible in the figure.
In some cases the anomalous dimension is even known exactly, e.g., for (partially) conserved currents.For the axialvector current the scale dependence, which we find from our nonperturbative calculations, is displayed in Fig. 10, where the renormalization condition (8) has been employed.Using instead Eq. ( 14) leads to a less satisfactory agreement with the expectations, although the conversion factor is exactly equal  to 1 in this case so that all sources of truncation errors disappear.While in Fig. 10 the deviation from 1 in the continuum limit remains below 0.002, it reaches about 0.007 in the same range of scales when Eq. ( 14) is applied.
Another case where truncation errors are absent in the anomalous dimension as well as in the conversion factor is given by the ratio Z S /Z P , see Fig. 11.
For the three-quark operators the conversion matrices are only known to one-loop accuracy, see Appendix G.As an example we consider the operator (B3) and compare in Fig. 12 Z(µ)/Z(µ 0 ) with the scale dependence obtained in threeloop perturbation theory.As above, we have used a secondorder polynomial in a 2 for the continuum extrapolation of Z(µ)/Z(µ 0 ).
As the figures show, there may be significant differences between the perturbative predictions for the scale dependence and the values of Z(µ)Z −1 (µ 0 ) even after an extrapolation to the continuum limit.When plotting Z RGI rather than these 12. Ratios of the renormalization factor of the three-quark operator (B3) evaluated in the RI ′ -SMOM scheme at different scales with the scale in the denominator fixed at µ 2 0 = 24 GeV 2 .The nonperturbative results at finite β have been computed in the RI ′ -SMOM scheme and a global chiral extrapolation has been performed.Finally, the values have been converted to the MS scheme using one-loop perturbation theory and extrapolated to the continuum limit (β = ∞).The curve shows the behavior calculated in three-loop perturbation theory in the MS scheme.ratios we should therefore expect to see deviations from a constant, not only due to the truncation of the perturbative series but also due to lattice artifacts.As an example we show in Fig. 13 Z RGI a2 computed with the help of the two-loop and the three-loop conversion factors from the RI ′ -MOM results.Using the locally chirally extrapolated data for the purpose of these plots enables us to reach rather large scales for the higher β values, where the lattice artifacts are particularly pronounced.
Indeed for most of the curves there is no interval in µ 2 where a plateau can be seen.There are several reasons for this behavior.At very small scales finite size effects might not be negligible and perturbation theory will probably break down completely.Moreover, the chiral extrapolation is less stable in this region.Beyond 2 or 3 GeV 2 truncation effects compete with lattice artifacts.Both effects tend to increase the Z RGI values (at least in the case at hand).The truncation effects are independent of β, decrease as a function of µ 2 and get smaller when higher orders in the perturbative expansion are taken into account.The effect of changing the loop order of the conversion factor can be seen by comparing the top and the bottom panels in Fig. 13.In both plots the scale dependence in the MS scheme has been computed from the fiveloop results for the anomalous dimension and the β function (see Appendix D).Using instead the four-loop approximation would only lead to hardly visible changes.
The discretization artifacts, on the other hand, depend only on a 2 µ 2 for fixed β, being proportional to a 2 µ 2 in a first approximation.They increase as a function of µ 2 , but decrease at a given µ 2 when β gets larger.The combination of both effects produces a structure in the data, which moves to  higher scales and eventually becomes a minimum as a gets smaller.Depending on the loop order there may appear "fake plateaus", e.g., for β = 3.34 in the two-loop case or for β = 3.55 when the three-loop conversion factor is employed.
A plot of Z RGI as a function of µ 2 for the three-quark operator (B3) is shown in Fig. 14.

VIII. TOWARDS FINAL RESULTS
In this section we present results for the renormalization and mixing coefficients obtained by two different methods.The first method, employed in our previous papers, e.g., Refs.[26][27][28], makes use of fits of the scale dependence.We call it the fit method in the following.In the second method we evaluate the renormalization factors at some suitable fixed scale and evolve them perturbatively to the desired reference scale.This method will be referred to as the fixedscale method.
In the fit method we try to exploit as much of the avail- able nonperturbative information as possible by performing a joint fit of the µ-dependence of the (chirally extrapolated) renormalization matrices Z(µ, a) MC for several β values.As explained above, the fit should describe (and hence suppress) the lattice artifacts vanishing like a power of a.
The choice of the fitting procedure is motivated by the following considerations.From the matrix of anomalous dimensions (evaluated for as many loops as possible) one can calculate a corresponding approximation of W (µ, µ 0 ) = Z(µ)Z −1 (µ 0 ), which should describe the µ-dependence for sufficiently large scales µ if there were no discretization effects.As our renormalization conditions respect the hypercubic symmetry of our lattices, the discretization artifacts must be H(4) invariant functions of the momenta defining the renormalization scheme.Invariant polynomials in the momentum components are, e.g., given by Hence an ansatz for the description of the lattice artifacts can be constructed from terms like see, e.g., Refs.[42][43][44].For the momenta p, q, . . .chosen in our renormalization conditions (see Eqs. ( 9), (10), and (51)) all these terms reduce to powers of a 2 µ 2 .
Adding such an ansatz for the effective description of the lattice artifacts, we arrive at the following fit function for the matrices Z(a, µ) MC : The fit parameters are the entries of the renormalization matrices Z(µ 0 , a) at the reference scale µ 0 and the entries of the three matrices A i parametrizing the lattice artifacts.Performing a joint fit for several β values we neglect a possible dependence of A i on the gauge coupling g 2 , which we only vary by 15%.
Representing W (µ, µ 0 ) as ∆Z(µ Therefore one has the equivalent options of using either Z(µ 0 , a) or Z RGI (a) as fit parameter.In the following we shall give our results for Z(µ 0 , a) with µ 2 0 = 4 GeV 2 , because these are more immediately useful, in particular in the presence of mixing.
The fluctuations at different scales are correlated.Hence we estimate the statistical uncertainty of the fit result Z(µ 0 , a) by the statistical error of the closest (in µ) data point.In most cases this procedure leads to errors that are still quite small.
The systematic uncertainties are much more important.In order to estimate them we perform a number of fits varying exactly one element of the analysis at a time, see Table III.We choose two values of the lower limit µ 1 of the fit range.For the parametrization of the discretization artifacts we either take the complete expression in Eq. ( 72) with n disc = 3 terms or we set A 3 = 0 corresponding to n disc = 2.The uncertainty due to the scale setting is taken into account by multiplying the values of 1/a 2 given in Table II by λ 2 scale = 1.03.This value contains the scale uncertainty of given in Ref. [21] and the largest error of our determination of t * 0 /a 2 , added in quadrature.Also Λ MS is varied within its uncertainty [21].In order to estimate the error caused by the truncation of the perturbative expansion of the conversion factors we reduce the number of loops n loops used in the calculation of the conversion factors by one, compared to the maximal value n max loops that is available.The size of the truncation error is taken to be one half of the resulting difference, because going, e.g., from two loops to three or more loops in the perturbative expansion is expected to lead to a smaller change than going from one loop to two loops, at least for sufficiently large scales.For the quark-antiquark operators we have n max loops = 3 in the RI ′ -MOM scheme, while in the RI ′ -SMOM scheme we have to be satisfied with n max loops = 2 for the tensor operators with derivatives.The symmetry properties of our operator multiplets can be found in Table IV in Appendix A, and numerical values of the coefficients in the perturbative expansion of the conversion factors are given in Table V in Appendix F. For the three-quark operators only one-loop results for the conversion factors are available so that in these cases n max loops = 1.Another source of systematic error is the chiral extrapolation, all the more as on the two finest lattices we have only two different masses at our disposal.We try to estimate the related uncertainty by comparing results obtained by means of the local and the global chiral extrapolation, cf.Sec.V.In the column of Table III labeled by χ these two choices are indicated by the letters g and l for the global and the local chiral extrapolation, respectively.It should however be noted that the upper limits of the fit ranges differ for data which have been locally or globally extrapolated.In both cases we include all available data points with µ ≥ µ 1 in the fit.In the case of the local chiral extrapolation the upper limit is therefore determined by the largest momenta considered.These have an approximately fixed value in lattice units corresponding to a 2 µ 2 ≈ 10 in the RI ′ -MOM scheme and a 2 µ 2 ≈ 5 in the RI ′ -SMOM scheme for the quark-antiquark operators and a 2 µ 2 ≈ 11 for the threequark operators.The global chiral extrapolation, on the other hand, requires the same scale in physical units for all β values.Thus the upper limit of the fit range is essentially determined by the largest scale available at the smallest β.In the case of the quark-antiquark operators, where this is β = 3.34, the fit range extends up to approximately µ 2 = 50 GeV 2 in the RI ′ -MOM scheme and µ 2 = 25 GeV 2 in the RI ′ -SMOM scheme.For the three-quark operators we do not have data at β = 3.34, and the fit range for the global chiral extrapolation extends up to approximately µ 2 = 75 GeV 2 .
As our central value we take the outcome of fit 1 with the statistical error determined as explained above.The six differences due to the discussed systematic uncertainties are added in quadrature to yield our final estimate of the systematic error.In most cases, the error due to the truncation of the perturbative expansion of the conversion factor is the largest contribution to the total systematic uncertainty.
The results obtained in this way, based on the RI ′ -MOM scheme, are collected in Tables VI-VIII, see Appendix H. Table VI (VII) contains results for operators with less than two derivatives obtained without (with) the perturbative subtraction of lattice artifacts.Here Z ′ V and Z ′ A have been determined with the help of the renormalization conditions ( 12) and ( 14), respectively, while for Z V and Z A the standard definition (8) has been used.In Table VIII we present the Z factors for operators with two derivatives, for which the perturbative subtraction of lattice artifacts is not yet available.
Ideally, Z factors determined with the help of the fit method should agree within errors, whether or not lattice artifacts have been perturbatively subtracted.A comparison of Tables VI and VII shows to which extent this expectation is fulfilled.Note that the fit method, which tries to suppress discretization effects as far as possible, is quite close in spirit to what is done in lattice perturbation theory, where (power-like) lattice artifacts are completely eliminated, cf.Sec.IV.
The fixed-scale method for the determination of the renormalization and mixing factors proceeds as follows.We first interpolate the chirally extrapolated data by cubic splines in ln(a 2 µ 2 ).Using this interpolation we can read off Z and its statistical error at some scale µ 2 , which we choose to be µ 2 = 10 GeV 2 , and evolve the result perturbatively to the desired reference scale µ 2 0 = 4 GeV 2 .Thus the statistical error stems from a higher scale than in the fit method and is therefore smaller.Again, the systematic uncertainties are usually considerably larger than the statistical error.The uncertainties due to the scale setting and the error in Λ MS are taken into account in the same way as in the fit method, cf. the entries for fits 4 and 5 in Table III.Also the error caused by the truncation of the perturbative expansion of the conversion factors is estimated through the same procedure as described above.The uncertainty due to the chiral extrapolation is again determined from the difference between the results obtained by means of the local and the global chiral extrapolation.These four systematic errors are added in quadrature to yield our final estimate of the total systematic uncertainty.The central value and its statistical error are taken from the interpolation of the globally chirally extrapolated data.
In Tables IX-XI (see Appendix H) we display the results coming from the fixed-scale method, based on the RI ′ -MOM scheme, separately for operators with less than two and with two derivatives obtained with and without the perturbative subtraction of lattice artifacts.In the fixed-scale approach we keep the lattice artifacts as they are generated in the simulations and try to get rid of them only when performing the continuum limit of physical (renormalized) quantities.Therefore, results obtained by this method with and without the perturbative subtraction of lattice artifacts need not coincide.By definition, the choice µ 2 = 10 GeV 2 is fixed and a variation of this value does not enter the systematics.
The corresponding results for the RI ′ -SMOM scheme are presented in Tables XII-XVII in Appendix H, again separately for the operators without derivatives (with and without perturbative subtraction of lattice artifacts) and for operators with derivatives.The results in Tables XII-XIV have been obtained with the help of the fit method, while results coming from the fixed-scale method can be found in Tables XV-XVII.In the cases where mixing occurs, i.e., for O v3 , O a2 , O h2a , O h 2b , and O h2c , we have given only the element Z 11 of the renormalization and mixing matrix, which can be compared with the RI ′ -MOM result.The full 2 × 2 mixing matrices can be found in ancillary files.In the case of the three-quark operators, where mixing matrices of size up to 4 × 4 appear, the results are presented in ancillary files only.
When evaluating renormalized quantities one can take into account the systematic uncertainties of the renormalization coefficients by error propagation and subsequent continuum extrapolation.However, in the presence of operator mixing, it might be more reasonable to use the various determinations of Z that go into the estimation of the systematic error, e.g., the results from fits 2-7 in the fit method, in order to compute the corresponding renormalized quantities and to use these numbers to estimate the systematic uncertainty of the renormalized quantity in the continuum limit.This is the procedure that we

FIG. 15. Ratios of different determinations of Z ′
A plotted against a 2 .For the black circles (red triangles) the fit (fixed-scale) method has been used, in the numerator within the RI ′ -MOM scheme and in the denominator within the RI ′ -SMOM scheme.The blue squares (green diamonds) represent results obtained in the RI ′ -MOM (RI ′ -SMOM) scheme with the fit method employed for the numerator and the fixed-scale method employed for the denominator.The systematic errors of the RI ′ -MOM numbers amount to 2 at most, while the RI ′ -SMOM numbers obtained with the fit (fixed-scale) method suffer from systematic uncertainties of up to 5 (2 ).
have applied in Refs.[26][27][28].Note that in Ref. [28] the renormalization and mixing coefficients were determined with the help of the fit method applied to a smaller set of locally chirally extrapolated data.In spite of these and a few further little differences the values used in Ref. [28] are consistent with our present results using two-loop conversion factors.An update of Ref. [28] using the new three-loop conversion factors can be found in Ref. [29].

IX. DISCUSSION
For a large set of quark-antiquark operators as well as for some three-quark operators we have presented renormalization factors computed with the help of various methods.For each of the quark-antiquark operators we end up with up to eight different results.When multiplied with the corresponding bare matrix elements, all of these renormalization factors should lead to the same continuum limit.This means that for a given operator, ratios Z (1) (µ, a)/Z (0) (µ, a) of two determinations Z (1) and Z (0) of the renormalization factor must tend to one as a → 0. However, one has to keep in mind that this statement holds only up to errors due to the truncation of the perturbative series entering the evaluation of Z (1) and Z (0) .The same considerations apply of course when one compares renormalization factors of a given operator computed using entirely different methods.
In the following we shall discuss a selection of examples of such ratios.In Figs. 15  .For the filled symbols the fixed-scale method has been used, in the numerator within the RI ′ -MOM scheme and in the denominator within the RI ′ -SMOM scheme.The open symbols represent results obtained in the RI ′ -SMOM scheme with the fit method employed for the numerator and the fixed-scale method employed for the denominator.The circles (triangles) show determinations using the threeloop (two-loop) conversion factor.For the results obtained with the three-loop conversion factors, the systematic uncertainty is up to 3 (7 ) in the case of the fixed-scale (fit) method.
different methods, plotted against a 2 .In all cases lattice artifacts have been subtracted perturbatively and the errors have been computed by error propagation of the statistical errors of the numerator and the denominator.These do not include the systematics discussed above.

Ratios of different determinations of Z ′
A are displayed in Fig. 15.As in this case the anomalous dimension vanishes identically and the conversion factor is exactly equal to one, there are no uncertainties due to truncations of perturbative expansions and all ratios considered here should tend to one in the continuum limit.The plot suggests that this is indeed the case within the statistical errors.
The next example (Fig. 16) is Z A .In this case the anomalous dimension is also known exactly, but the conversion factor is only known to three loops.Therefore truncation errors are to be expected.Indeed, the results obtained with the threeloop conversion factor differ significantly from those obtained with the two-loop conversion factor and extrapolate to a value closer to one.
A similar behavior is found for operators with one derivative.For these operators also the anomalous dimension is only known to a finite order in perturbation theory.However, it turns out that the perturbative expansion of the anomalous dimension converges much faster than the expansion of the conversion factor.As an example we consider Z v 2b in Fig. 17, where the five-loop expression for the anomalous dimension is used throughout, but the number of loops taken into account in the conversion factor is varied.The effect of this variation is clearly visible and goes into the desired direction.
Finally, we show in Fig. 18  .For the filled symbols the fit method has been used, in the numerator within the RI ′ -MOM scheme and in the denominator within the RI ′ -SMOM scheme.The open symbols represent results obtained in the RI ′ -MOM scheme with the fit method employed for the numerator and the fixed-scale method employed for the denominator.The circles (triangles) show determinations using the three-loop (twoloop) conversion factor.For the results obtained with the three-loop conversion factors, the systematic uncertainty amounts to about 1% (2-3%) in the case of the RI ′ -MOM (RI ′ -SMOM) scheme.
tions of Z S , with the RI ′ -MOM (RI ′ -SMOM) scheme employed in the numerator (denominator).For a → 0 the value of Z S determined with the help of the RI ′ -SMOM scheme is about 5% larger than the number from the RI ′ -MOM scheme.Unusually large differences between the two methods for the scalar density have already been reported in Ref. [45] at two relatively large values of the lattice spacing, using a different lattice action.In that case the matching between the RI ′ -SMOM and the MS scheme could only be carried out at twoloop order.In Fig. 18 we only display the statistical errors, however, also when including the systematic uncertainties, the result is unsatisfactory.This becomes even more apparent for the outcome of the fixed-scale method, which carries smaller statistical and systematic errors.In both cases (fit and fixed-scale method) the results obtained with the three-loop conversion factor are closer to unity than those obtained with the two-loop conversion factor.In comparison to the other operators, the perturbative coefficients for the matching between the MS and the RI ′ -SMOM schemes are quite small, see Table V.Therefore, we suspect that our estimate of the perturbative uncertainty as half of the difference between the results obtained from two-and three-loop matching may be an underestimate for this particular operator.Unfortunately, we cannot carry out a similar comparison between the RI ′ -SMOM and the RI ′ -MOM schemes for Z P , which shares its perturbative matching and its anomalous dimension with Z S , since the RI ′ -MOM scheme is unsuitable in this case.
Within CLS some renormalization factors have also been computed by other groups, utilizing the RI ′ -MOM scheme [10]  .For the filled (open) symbols the fit (fixed-scale) method has been used, in the numerator within the RI ′ -MOM scheme and in the denominator within the RI ′ -SMOM scheme.The circles (triangles) show determinations using the three-loop (two-loop) conversion factor.For the results obtained with the three-loop conversion factors, the systematic uncertainty amounts to 2% (3-4%) in the case of the RI ′ -MOM scheme and the fixed-scale (fit) method, it is about 0.4% (0.5-1.5%) in the case of the RI ′ -SMOM scheme and the fixed-scale (fit) method.14,15].In Figs.19 -22 we show for a few examples ratios of such alternative results divided by our numbers.For this comparison we employ our results with the perturbative subtraction of lattice artifacts and add our statistical and systematic errors in quadrature, since we cannot separate these two sources of uncertainties in the case of the alternative determinations.
We begin with Z S /Z P and Z P /(Z S Z ′ A ).In these cases, our results are free of perturbative ambiguities, because the anomalous dimensions as well as the conversion factors cancel between numerator and denominator.Nevertheless, systematic errors can be rather large, as, e.g., Table XII shows.We denote the renormalization factor of the axialvector current by Z ′ A , because for our numbers we use the data obtained with the help of the renormalization condition (14).The ratio Z S /Z P has been studied by the ALPHA collaboration in Ref. [15].We divide the values given in the column "WI(1468)" of Table 6 in Ref. [15] by our numbers determined with the help of the fixed-scale and the fit method.The results are shown in Fig. 19 as black circles and red triangles, respectively, plotted against a 2 .Indeed, these ratios of independent determinations appear to approach unity in the continuum limit.However, the curvature that we see when plotting the data as a function of a 2 indicates that within our range of lattice spacings the discretization effects are not purely O(a 2 ).It needs to be seen which of the three methods (Ref.[15], fit method, fixed-scale method) will result in the most convincing continuum limit extrapolation, once Z S /Z P is multiplied by corresponding ratios of matrix elements, computed in lattice simulations.[15] divided by our results obtained in the RI ′ -SMOM scheme with the fixed-scale (fit) method.
In Fig. 20 we divide the combination Z P /(Z S Z A ) of the ALPHA collaboration given in Table 4 (trajectory LCP-1) of Ref. [14] by our results on Z P /(Z S Z ′ A ) obtained in the RI ′ -SMOM scheme with the fixed-scale method (black circles).The data appear to overshoot the continuum limit expectation by about 2%.Repeating this exercise for the Z P /(Z S Z A ) combination, determined by the RQCD collaboration [46] from a fit to axial Ward identity quark masses [5] (red triangles), some curvature of the data becomes apparent.Still, the figure does not contradict the expectation that the value of one will be reached in the continuum limit.
In Fig. 21 we plot ratios of the renormalization factor Z A determined by the ALPHA collaboration [12], divided by our results on Z ′ A , obtained with the fixed-scale and the fit methods.Both ratios tend to unity in the continuum limit and again this approach is not linear in a 2 within the range of lattice spacings covered.We can also compare the ratio Z A /Z P of Ref. [13] with our results on Z ′ A /Z P obtained from the RI ′ -SMOM scheme.Using either the fit or the fixed-scale method, in the continuum our results appear to be smaller by 2%.In the latter case this deviation exceeds the combined statistical and systematic errors of our numbers by a factor of about three.Given that the axialvector renormalization constants agree between different determinations, as does Z S /Z P (see Figs. 21  and 19), this tension is consistent with the observation made in Fig. 18 that Z S is larger using the RI ′ -SMOM scheme than using the RI ′ -MOM scheme.Therefore, we cannot exclude the possibility that our Z P is somewhat overestimated, e.g., due to perturbation theory uncertainties.
As our last example we consider Z v 2b , the renormalization factor of an operator containing a covariant derivative.In Fig. 22 our results are compared with the numbers obtained in Ref. [10].These were determined in the RI ′ -MOM scheme at β = 3.40, 3.46 and 3.55 using an approach similar to our fit method, though some details differ.Therefore we compare  4 (trajectory LCP-1) in Ref. [14] divided by our results obtained in the RI ′ -SMOM scheme with the fixed-scale method.The red triangles show data of the RQCD collaboration [46] for ZP /(ZSZA) determined from axial Ward identity masses [5] divided by our results obtained in the RI ′ -SMOM scheme with the fixed-scale method.

FIG. 21. Ratios of different determinations of Z ′
A plotted against a 2 .The black circles (red triangles) show the values given in the column Z l A,sub of Table 7 in Ref. [12] divided by our results obtained in the RI ′ -SMOM scheme with the fixed-scale (fit) method.them with our fit results obtained in the RI ′ -MOM scheme.Within the errors the ratios are compatible with unity for all three β values, although the numbers given in Ref. [10] seem to lie consistently above ours.Similar observations apply for the other operators with one derivative studied in Ref. [10].
In view of the multitude of up to eight determinations of the renormalization factors for a single operator given in this paper, the question arises, which values should be used in a particular situation.To answer this question a few criteria can be given.The decision in favor of the RI ′ -MOM or RI ′ -SMOM intermediate scheme is obvious when nonforward ma- trix elements of operators are required which mix with totalderivative operators.In such cases the RI ′ -SMOM scheme is mandatory.It is also strongly preferred for the pseudoscalar density, because the pion pole makes the chiral extrapolation problematic in the RI ′ -MOM scheme.For all operators that are unaffected by the pion pole and for which the kinematics of the RI ′ -MOM scheme is admissible, the choice of the intermediate scheme is purely a matter of taste, with the possible exception of the scalar density, where the RI ′ -SMOM scheme may be more favorable because of the smaller coefficients in the perturbative expansion (68) of the conversion factor (see Table V).
Concerning the distinction between the fit method and the fixed-scale method one can say that the fit method tries to suppress the power-like lattice artifacts as far as possible.The fixed-scale method, on the other hand, does not care about lattice artifacts so that they will be taken into account only in the continuum extrapolation of the physical observables.If the data for the observables allow to perform such an extrapolation reliably, the fixed-scale method could be preferable.However, if data are available for a few lattice spacings only (maybe only for a single value) and a decent continuum extrapolation is impossible, renormalization factors determined with the help of the fit method might lead to estimated results which are closer to the continuum value.Similarly, the perturbative subtraction of lattice artifacts may be more useful when a careful continuum limit is still out of reach, while it does not help that much when the power-like lattice artifacts are taken care of by the continuum extrapolation anyway.

X. SUMMARY
Using the Rome-Southampton method and variants thereof, we have computed renormalization and mixing factors nonperturbatively within the CLS setup, i.e., for n f = 2 + 1 flavors of nonperturbatively improved Wilson quarks and the Lüscher-Weisz action with tree-level coefficients for the gluons.We have developed an approach that allows us to include in our analysis also ensembles with open boundary conditions in the time direction.Quark-antiquark operators as well as three-quark operators have been considered with the help of a variety of methods.The fit method has already been employed previously in our papers [26,28].Comparing the results obtained with different methods allows us to check the consistency of the employed procedures in the continuum limit.For applications we recommend the numbers from the fixed-scale method given in Tables XVI and XVII in Appendix H.The corresponding values from the fit method are collected in Tables XIII and XIV.The 2 × 2 mixing matrices of quark-antiquark operators and the results for the threequark operators can be found in ancillary files.
Regarding the renormalization factor Z S , we see inconsistencies on the several percent level that at present we cannot explain.Similar problems have been reported in Ref. [45].Since the ratio Z S /Z P appears to be consistent with other determinations and the discrepancy reduces if the matching to the MS scheme is carried out at a higher perturbative order, we suspect that for the scalar and pseudoscalar densities, there may be an unusually large perturbation theory uncertainty of about 2-3% that we underestimated.By implication, the results on Z ′ A /Z P should be considered with caution.In many cases, with our accuracy and within the range of lattice spacings covered, the approach to the continuum limit of ratios between determinations using different methods is not a linear function of a 2 .Therefore, it is advisable to employ more than one set of renormalization factors in continuum extrapolations of nonperturbative matrix elements.
At present, on the two finest lattices there are only two different masses each available.This fact limits the accuracy of the required chiral extrapolations.Generating additional ensembles with smaller masses is desirable.
Another limitation of the accuracy of our results arises from the necessity to use continuum perturbation theory.For many operators the conversion or matching factors needed to obtain results in the MS scheme are now known to three loops, but the convergence of these perturbative expansions is generally not too good, in contrast to the anomalous dimensions, which control the scale dependence.Therefore further efforts in the perturbative evaluation of matching factors would be helpful, especially for the (pseudo)scalar operator, but also for threequark operators, for which the matching factors are presently only known to one-loop accuracy.The authors thank Oleg Veretin and John Gracey for their valuable calculations in continuum perturbation theory.
We used a modified version of the CHROMA [47] software package along with improved inverters [48][49][50][51].The configurations were generated as part of the CLS effort [3,5] using OPENQCD [48,52].We thank all our CLS colleagues for the joint generation of the gauge ensembles.Additional ensembles with degenerate light and strange quark masses were generated with OPENQCD by members of the Mainz group on HPC Clusters of IKP Mainz as well as by RQCD on the QPACE systems of the SFB/TRR 55 using the BQCD code [53].Correlation functions were computed on the QPACE 3 machine of the SFB/TRR 55.
Computer time on the DFG-funded Ara cluster at the Friedrich-Schiller-University Jena is acknowledged.A. Sternbeck acknowledges support by the BMBF under Grant No. 05P15SJFAA (FAIR-APPA-SPARC) and by the DFG Research Training Group GRK1523.
The authors gratefully acknowledge the Gauss Centre for Supercomputing (GCS) for providing computing time for GCS Large-Scale Projects on the GCS share of the Su-perMUC system at Leibniz Supercomputing Centre (LRZ, https://www.lrz.de) and JUWELS [54] at Jülich Supercomputing Centre (JSC, http://www.fz-juelich.de/ias/jsc/).GCS is the alliance of the three national supercomputing centres HLRS (Universität Stuttgart), JSC (Forschungszentrum Jülich), and LRZ (Bayerische Akademie der Wissenschaften), funded by the German Federal Ministry of Education and Research (BMBF) and the German State Ministries for Research of Baden-Württemberg (MWK), Bayern (StMWFK) and Nordrhein-Westfalen (MIWF).
Appendix A: Quark-antiquark operator multiplets In this Appendix we list the multiplets of quark-antiquark operators that we consider.Flavor indices are omitted for simplicity.Color and spinor indices are suppressed.As usual, {• • • } will denote the symmetrization of all enclosed indices.While a universal factor multiplying all members of a multiplet is irrelevant for the renormalization, it is important to ensure that the ratios of the individual normalization factors are such that the operators transform under H(4) according to a unitary representation.
We start with the operators without derivatives.Their mass dimension equals three and they do not mix with other operators (operator multiplets).The scalar density and the pseudoscalar density form multiplets of dimension one.In the cases of the local vector current and the local axialvector current we have multiplets of dimension four (µ = 1, 2, 3, 4).The components of the tensor density (or tensor current) with 1 ≤ µ < ν ≤ 4 make up a six-dimensional multiplet.
The operators with a single covariant derivative have mass dimension four.We consider the following multiplets, which in the flavor-nonsinglet case do not mix with any other operators.The names we give to the operators and to the corresponding renormalization factors (matrices) are motivated by the nomenclature for the moments of the parton distribution functions of the nucleon.
The v 2a operators where 1 ≤ µ < ν ≤ 4, form a multiplet of dimension six, while the v 2b operators span a three-dimensional multiplet.
Analogously we define the r 2a operators with 1 ≤ µ < ν ≤ 4, and the r 2b operators With the help of the abbreviation we can write the eight h 1a operators as

(A11)
The h 1b operators form an eight-dimensional multiplet as well.
The operators with two derivatives are of mass dimension five.In the flavor-nonsinglet case at leading twist one can choose them such that they do not mix with operators of lower mass dimension.However, when nonforward matrix elements are needed, as they appear, e.g., in the RI ′ -SMOM scheme, we cannot avoid the mixing with the so-called total-derivative operators, which occurs already in the continuum.Therefore we consider in the case of the v 3 operators the two fourdimensional multiplets with {µ, ν, λ, ρ} = {1, 2, 3, 4}.They transform identically under H(4).Note that in the continuum the operator O (ρ) 1 can be written as 2 is the second derivative of the vector current: On the lattice, this relation is violated because of lattice artifacts in the derivatives.Analogous remarks hold also for the remaining operators with two derivatives.These comprise the a 2 operators = ψ(x) again with {µ, ν, λ, ρ} = {1, 2, 3, 4}, as well as the h 2a , h 2b , and h 2c operators.
In order to define the latter we use the abbreviation The first multiplet of the h 2a operators is then given by Replacing O T − by O T + yields the corresponding second multiplet O 2 with the identical transformation behavior under H(4).In the case of the h 2b operators we have the first multiplet Again, the second multiplet O 2 is obtained by replacing O T − with O T + .Finally, as the first multiplet of the h 2c operators we take The corresponding multiplet O 2 is constructed as in the two previous cases.

+1
Or 2b τ (3) 4 In the RI ′ -MOM case the respective second multiplets O 2 do not contribute because their forward matrix elements vanish.
In Table IV we give the H(4) irreducible representations in the notation of Ref. [55] and the charge conjugation parity for our operator multiplets.
We summarize in this Appendix our results for the term F (a 2 µ 2 ) in Eq. ( 52).This term is calculated numerically in lattice perturbation theory for the quark propagator and all quark-antiquark operators with less than two derivatives for a set of renormalization scales a 2 µ 2 .To be consistent with our nonperturbative calculation, we use the Landau gauge, apply exactly the RI ′ -(S)MOM renormalization conditions described in section II B and use the lattice expressions for the propagators and vertices corresponding to our action.In particular, the improved gauge action [56] leads to a complicated gluon propagator consisting of many terms, which slows down the calculation considerably.
The analytical part of the calculation is performed in FORM (v4.2.1) [57].In FORM we implement the RI ′ -(S)MOM renormalization conditions and insert the one-loop expres-sions for the respective correlation functions.The corresponding diagrams are shown in Fig. 23.For the quark self-energy there are two diagrams, a sunset and a tadpole diagram, and for all operator insertions a vertex diagram.In addition there is a tadpole and a leftand a right-sail diagram for operators with one derivative.
The propagators and the quark-gluon vertex on the lattice can be found in standard textbooks or in the review [58].This review contains many useful expressions, for example, the gluon propagator and the three vertices for the operator insertion ψ(x)γ µ → D ν ψ(x) at zero momentum transfer, which enter the calculation of the RI ′ -MOM renormalization constants Z v2a and Z v 2b .
For the RI ′ -SMOM scheme the operator vertices at finite momentum transfer are required.A straightforward calculation, e.g., for the operator ψ(x)γ µ ↔ D ν ψ(x) yields contributions of order g 0 , g 1 and g 2 : For brevity, the factor (2π) 4 and the momentum conserving δ function have been omitted.The incoming and outgoing (a) 23. One-loop diagrams for the quark self-energy (a,b) and quark bilinears (c,d,e,f) with momentum transfer at the operator insertion (solid box).For operators without derivative the last row of diagrams is absent.) Numerical results for F (a 2 µ 2 ) versus a 2 µ 2 for operators with (bottom) and without (top) derivative in the RI ′ -MOM scheme (left panels) and in the RI ′ -SMOM scheme (right panels), displayed as data points connected by lines.The top panels also show F (a 2 µ 2 ) for the quark propagator where the quark momentum direction matches that of the respective RI ′ -(S)MOM condition.Analytical results for F (a 2 µ 2 ) valid up to O(a 2 ) are shown as solid lines.The deviations from zero of the data points at small a 2 µ 2 are numerical artifacts and can be ignored.
quark momenta are denoted by p and p ′ , respectively, and Q = p ′ − p is the quark momentum transfer at the vertex.The matrices T a are the standard SU(3) generators.For operators without derivatives, there is only V (0) and this equals the gamma matrix structure of the respective operator insertion, because the sine and cosine functions in Eq. (C1) are due to the derivative in the operator.Taking the traces in the renormalization conditions, FORM generates very long expressions for the one-loop part of the renormalization constants.These expressions have to be integrated over the internal gluon loop momentum (see Fig. 23).This is performed numerically using a separate C++ program and the library cubature [59].The external momenta are kept fixed, i.e., for each setup of external quark momenta a separate numerical integration is performed.The quark momenta are chosen to exactly match those for the nonperturbative calculation.
To speed up the integration, we use FORM's output optimization feature (Format O3) to generate compact integration kernels as input for our C++ program.Compared to the other FORM operations, this final output optimization is CPU-time intensive and uses a Monte Carlo Tree Search to find an optimal Horner scheme for a kernel.In this way, the final expressions for the integration kernels become compact multivariate polynomials whose numerical evaluation (for each integration step) is less CPU-time intensive than that of the original kernel.
After integration the one-loop contribution to Z for each operator can be written as (cf.Eq. ( 52)) where the first two terms are known from lattice perturbation theory in the limit a → 0 and F (0) = 0. Subtracting these terms from the numerical results for z(a 2 µ 2 ) gives F (a 2 µ 2 ).
Note that F at fixed a 2 µ 2 would also depend on the directions of the external momenta if they are varied.For our RI ′ -(S)MOM renormalization conditions, however, the momentum directions are fixed, cf.Eqs. ( 9) and (10).Hence there is only a single curve F (a 2 µ 2 ) for each operator.Fig. 24 shows the numerical results for F (a 2 µ 2 ), separately for the RI ′ -MOM and RI ′ -SMOM schemes as well as for operators without derivatives and with one derivative.In the case of the RI ′ -MOM scheme, analytical expressions for F (a 2 µ 2 ) valid up to O(a 2 ) can be derived from results in the literature [35,60].They are included in Fig. 24 for comparison (solid lines) and agree well with our numerical results for a 2 µ 2 ≲ 0.4.Beyond this point deviations are clearly seen and corrections of higher order in a seem to become important.

Appendix D: Renormalization group functions for quark-antiquark operators
For the reader's convenience we collect in this Appendix perturbative results for the QCD β function and the anomalous dimension of the quark-antiquark operators considered in this paper.
In the MS scheme the coefficients of the β function, as defined in Eqs. ( 58) and ( 59), are given by [61][62][63]  Here and in the following, n f denotes the number of flavors and ζ n is the value of Riemann's ζ-function at n.
We now turn to the coefficients of the anomalous dimension in the MS scheme employing the conventions given in Eqs. ( 61) and ( 62).An anticommuting γ 5 is assumed.
Of course, the vector current and the (nonsinglet) axialvector current have vanishing anomalous dimension.For T µν we have [64] For S and P one finds [65,66] The five-loop coefficient has been calculated as well [67]: Note that the coefficient of n f in γ 2 was incorrectly given in Ref. [23] due to a misinterpretation of the results of Ref. [73].The ensuing difference in the anomalous dimension is however small, because the coefficient changes just from 392.948 to 399.275.The four-loop coefficient γ 3 can be found in Ref. [74]: The three-loop anomalous dimension of the h 1a and h 1b operators can be found in Refs.[75,76]: Finally, we can take the anomalous dimension of the quark field (see Eq. ( 63)) in Landau gauge to four loops from Ref. [77]: In the case of the RI ′ -SMOM scheme we have to take into account that operators with two derivatives mix with totalderivative operators.This concerns the v 3 and a 2 operators as well as the h 2a , h 2b , and h 2c operators.The anomalous dimension becomes a 2 × 2 matrix, whose 1-1 element coincides with the anomalous dimension given above.For the v 3 and a 2 operators the only other nonzero entry is the 1-2 element.This matrix element can be calculated from the results given in Refs.[76,78].Transforming these results to our operator basis one finds In our computations we use also the three-loop coefficient (γ 2 ) 12 , which can be evaluated with the help of unpublished results of John Gracey [79].Recently, the three-loop coefficients missing in Ref. [78] have been calculated numerically [8].Utilizing these numbers one finds (γ 2 ) 12 = −417.165+ 64.0972 n f + 1.09319 n 2 f .(D36) For the h 2a , h 2b , and h 2c operators the anomalous dimension matrix has three nonvanishing entries.The 1-1 element is identical to the anomalous dimension of these operators given above, while the 2-2 element coincides with the anomalous dimension of the tensor density T µν .Transforming the results given in Ref. [27] to our operator basis one obtains for the last nonvanishing matrix element In this Appendix we present the anomalous dimensions of the three-quark operators that we consider.For the operators without derivatives three-loop results have been obtained in Ref. [34].
For the multiplets (B1) and (B2) transforming according to the H(4) representation τ For the operators with one derivative, only one-loop anomalous dimensions are available.For the multiplet (B6) we have The 4 × 4 anomalous dimension matrix of the octet multiplets (B7) -( B10) is diagonal in one-loop order: In the case of the decuplet multiplets (B11) -(B13) we get a 3 × 3 anomalous dimension matrix: Appendix F: Conversion factors for quark-antiquark operators In this Appendix we compile results for the expansion coefficients of the conversion matrices (factors) of our operators.The conversion matrices C are defined in Eq. (17) and their perturbative expansion is given in Eq. (68).
We start with the RI ′ -MOM scheme.From Ref. [77]

(F9)
Perturbative expressions for the RI ′ -SMOM vertex functions are available in analytic form for up to two loops [27,78,81,82].As they are rather complicated, also the corresponding results for the coefficients of the conversion factors as derived from Eq. ( 17) are quite lengthy.Therefore we give them only in numerical form for n f = 3, see Table V. Threeloop results for the operators without derivatives have recently been obtained in Ref. [6] and for some operators with derivatives in Ref. [8].In these papers the required integrals have been evaluated with the help of numerical integration.The corresponding results for c 3 computed from the given central values are also included in Table V.Their numerical uncertainty is calculated by means of error propagation from the errors of the individual contributions given in Refs.[6,8].For our final renormalization factors this uncertainty is completely irrelevant.In the case of the scalar density, there is even an analytic expression for c 3 [7], which agrees perfectly with the numerical result.
Note that here the mixing with total-derivative operators has been ignored.Hence the numbers in Table V for the operators with two derivatives are sufficient only when forward matrix elements are considered.In order to facilitate the comparison between the RI ′ -SMOM and the RI ′ -MOM scheme we include also the numerical values of the coefficients in the RI ′ -MOM scheme.While for the scalar and the pseudoscalar densities the coefficients in the RI ′ -SMOM scheme are substantially smaller than in the RI ′ -MOM scheme, they are quite similar in the other cases.
When nonforward matrix elements are to be investigated, we have to take into account that operators with two deriva- In this Appendix we present the conversion matrices (factors) of the three-quark operators that we consider, as obtained in Ref. [33] to one-loop accuracy.Here we need the trigamma function with the special values ψ 1 (1/3) ≈ 10.0956 and ψ 1 (1/4) ≈ 17.1973.

ADDENDUM
After the publication of our article we noticed that some results from the literature used were incorrect.This affects in particular the three-loop conversion coefficient for the tensor density between the RI ′ -MOM and the MS schemes and the subtraction of lattice artifacts at order g 2 in the coupling for our operators of dimension three and four.The resulting changes are below 1 and smaller than the systematic uncertainties given, in most cases even considerably smaller.This will not lead to relevant differences in phenomenological results.Nevertheless, we prefer to correct the affected tables accordingly.We also include the four-and five-loop anomalous dimensions of the tensor density and the quark field, respectively, in the determination.The details are as follows.
In the perturbative subtraction of lattice artifacts we need the gluon propagator for the Lüscher-Weisz action with treelevel coefficients.This has been taken from Eqs. (11.26) - (11.28) in Ref. [58].Unfortunately, there is an error in Eq. (11.28) as can be seen from the original literature [56,83].Moreover, the conversion factor (up to three loops) for the tensor density T µν in the RI ′ -MOM scheme was extracted from results given in Ref. [80].In Ref. [84] it is pointed out that the three-loop contribution is incorrect.The coefficient c 3 should read For n f = 3 one finds c 3 = −1194.56to be compared with the value c 3 = −1207.96given in Table V .We also realized that the MS anomalous dimension of T µν is actually known to four-loop accuracy [85], while in our paper we used only the three-loop anomalous dimension.The four-loop coefficient reads (Z.4) We have recalculated the lattice artifacts in one-loop lattice perturbation theory for the quark-antiquark operators with less than two derivatives using the correct gluon propagator and implemented the above conversion factor for T µν in the RI ′ -MOM scheme as well as the four-loop anomalous dimension of T µν and the five-loop anomalous dimension of the quark field.
With all these corrections in place, our analysis yields the Tables XVIII -XXV, where we include also the numbers that have not changed.
Tables XVIII, XIX, XX and XXI are the corrected versions of Tables VI, VII, IX and X, respectively.They differ from the latter tables due to the corrected subtraction, the updated anomalous dimensions of T µν and the quark field as well as due to the corrected value of the three-loop coefficient in the conversion factor for T µν .
Tables XXII, XXIII, XXIV and XXV are the corrected versions of Tables XII, XIII, XV and XVI, respectively.They differ from the latter tables due to the corrected subtraction as well as due to the updated anomalous dimensions of T µν and the quark field.The above-mentioned corrections and updates are contained in the Erratum to our paper [87].We now present further improvements resulting from the fact that for some operators the conversion factors leading from the RI ′ -MOM scheme to the MS scheme are meanwhile available to us in four-loop order, while in Ref. [87] only the three-loop conversion fac-tors were used.We previously overlooked that we could have extracted the four-loop expression for the conversion factor of the quark wave function renormalization C q from Ref. [88], and recently the four-loop coefficient in the conversion factor for T µν has been published in Ref. [84], while for S, V µ and A µ it can be extracted from Ref. [89].
One finds in the case of the quark wave function renormalization The numerical values are 188150.13for S, −14625.21for V µ and A µ , and −37788.57for T µν .
In Tables XXVI -XXIX we present the results obtained using the four-loop conversion factors for S, V µ , A µ , T µν and for the quark wave function renormalization.Correspondingly, the uncertainty due to the truncation of the expansion of the conversion factors is estimated by comparing with the three-loop expressions.The other elements of the evaluation are the same as in Tables XVIII -XXI.
The changes caused by the corrected subtraction, the updated anomalous dimensions and the corrected value of the three-loop coefficient in the conversion factor for T µν look rather small.However, replacing for Z q , Z S , Z V , Z A and Z T the three-loop (two-loop) RI ′ -MOM conversion factors by the four-loop (three-loop) expressions leads to larger differences.

FIG. 7 .
FIG.7.Ratios of the renormalization factor of the tensor density evaluated in the RI ′ -SMOM scheme at different scales with the scale in the denominator fixed at µ 2 0 = 24 GeV 2 .The nonperturbative results at finite β have been computed in the RI ′ -SMOM scheme, lattice artifacts have been subtracted perturbatively, and a global chiral extrapolation has been performed.Finally, the values have been converted to the MS scheme using one-loop, two-loop and three-loop (from top to bottom) perturbation theory and extrapolated to the continuum limit (β = ∞).The curves show the behavior calculated in three-loop perturbation theory in the MS scheme.

FIG. 8 .FIG. 9 .
FIG.8.The same as the lowest panel in Fig.7, but without the perturbative subtraction of lattice artifacts.

FIG. 10 .
FIG.10.The same as the lowest panel in Fig.7, but for the axialvector with the renormalization condition(8).

FIG. 11 .
FIG.11.The same as the lowest panel in Fig.7, but for the ratio ZS/ZP .
FIG.12.Ratios of the renormalization factor of the three-quark operator (B3) evaluated in the RI ′ -SMOM scheme at different scales with the scale in the denominator fixed at µ 2 0 = 24 GeV 2 .The nonperturbative results at finite β have been computed in the RI ′ -SMOM scheme and a global chiral extrapolation has been performed.Finally, the values have been converted to the MS scheme using one-loop perturbation theory and extrapolated to the continuum limit (β = ∞).The curve shows the behavior calculated in three-loop perturbation theory in the MS scheme.

85 FIG. 13 .
FIG.13.Z RGI for Oa 2 computed from the locally chirally extrapolated RI ′ -MOM results with the help of the two-loop (top panel) and the three-loop (bottom panel) conversion factor.In both cases the five-loop anomalous dimension has been used.

85 FIG. 14 .
FIG.14.Z RGI for the three-quark operator (B3) computed from the locally chirally extrapolated RI ′ -SMOM results with the help of the one-loop conversion factor.For the anomalous dimension the threeloop approximation has been used.

ACKNOWLEDGMENTS
This work was funded in part by the Deutsche Forschungsgemeinschaft (collaborative research centre SFB/TRR-55) and by the European Union's Horizon 2020 Research and Innovation programme under the Marie Skłodowska-Curie grant agreement no.813942 (ITN EuroPLEx).

TABLE I .
List of ensembles.The inverse gauge coupling β determines the lattice spacing, while the spatial and temporal extents fix the lattice geometry N 3 s × Nt.Boundary conditions in the time direction are either periodic (p) or open (o).The hopping parameter κ determines the corresponding quark mass; the resulting approximate pion mass mπ is given in units of MeV, followed by the spatial lattice size in pion mass units.
FIG.1.Components of the quark propagator on ensemble J500 summed over space as a function of time (in lattice units).The magnitude of the momentum is given by µ ≈ 2 GeV and a −1 ≈ 5 GeV.The volume for the momentum source is limited by the dotted grey lines.The volume for the final summation is indicated by the blue lines.

TABLE III .
Choices for the fits in the case n max loops = 3.Further explanations are given in the text.
-18 we show for a few operators ratios of renormalization factors determined in this paper by or Schrödinger functional techniques[11, 12, . The black circles (red triangles) show the values of ZS/ZP given in the column "WI(1468)" of Table6in Ref.
FIG. 20.Ratios of different determinations of ZP /(ZSZ ′A ) plotted against a 2 .The black circles represent the numbers for ZP /(ZSZA) in Table [10]22.Ratios of different determinations of Zv 2b plotted against a 2 .The values given in Table4of Ref.[10]have been divided by our results obtained in the RI ′ -MOM scheme with the fit method.

TABLE IV .
Operator multiplets and their transformation behavior under H(4).The charge conjugation parity is denoted by C.
[71,72]a , v 2b , r 2a , and r 2b operators have the same anomalous dimension.From Ref.[68]we obtain For the v 3 and a 2 operators we extract from Refs.[71,72] [80] the help of Ref.[80]we find for T µν (for a correction, see Eq. (Z.2) of the Addendum)

TABLE V .
Coefficients of conversion factors for n f = 3.In the cases where loop integrals have been evaluated numerically estimates of the resulting uncertainties have been included.For T in the MOM scheme c3 should read −1194.56instead of −1207.96,see the Addendum.This concerns the v 3 and a 2 operators as well as the h 2a , h 2b , and h 2c operators.One finds for the v 3 and a 2 operators the matrices tives mix with total-derivative operators.In this case we must use the RI ′ -SMOM scheme, and the conversion factors become conversion matrices of size 2 × 2 whose 1-1 entry coincides with what one gets from TableV.Appendix G: Conversion factors for three-quark operators

TABLE VIII .
Results from fits for operators with two derivatives based on the RI ′ -MOM scheme without the perturbative subtraction of lattice artifacts.The chiral extrapolation has been performed globally.The first number in parentheses gives the statistical error, while the second number is an estimate of the systematic uncertainty.All values refer to the MS scheme at the scale µ 2 0 = 4 GeV 2 .

TABLE X .
Results for operators with less than two derivatives based on the RI ′ -MOM scheme, obtained by means of the fixed-scale method with the perturbative subtraction of lattice artifacts.The chiral extrapolation has been performed globally.The first number in parentheses gives the statistical error, while the second number is an estimate of the systematic uncertainty.All values refer to the MS scheme at the scale µ 2 0 = 4 GeV 2 .See Table XXI of the Addendum for updated results.

TABLE XI .
Results for operators with two derivatives based on the RI ′ -MOM scheme, obtained by means of the fixed-scale method without the perturbative subtraction of lattice artifacts.The chiral extrapolation has been performed globally.The first number in parentheses gives the statistical error, while the second number is an estimate of the systematic uncertainty.All values refer to the MS scheme at the scale

TABLE XIII .
Results from fits for operators with less than two derivatives based on the RI ′ -SMOM scheme with the perturbative subtraction of lattice artifacts.The chiral extrapolation has been performed globally.The first number in parentheses gives the statistical error, while the second number is an estimate of the systematic uncertainty.All values refer to the MS scheme at the scale µ 2 0 = 4 GeV 2 .See Table XXIII of the Addendum for updated results.

TABLE XIV .
Results from fits for operators with two derivatives based on the RI ′ -SMOM scheme without the perturbative subtraction of lattice artifacts.The chiral extrapolation has been performed globally.The first number in parentheses gives the statistical error, while the second number is an estimate of the systematic uncertainty.All values refer to the MS scheme at the scale µ 2 0 = 4 GeV 2 .

TABLE XV .
Results for operators with less than two derivatives based on the RI ′ -SMOM scheme, obtained by means of the fixed-scale method without the perturbative subtraction of lattice artifacts.The chiral extrapolation has been performed globally.The first number in parentheses gives the statistical error, while the second number is an estimate of the systematic uncertainty.All values refer to the MS scheme at the scale µ 2 0 = 4 GeV 2 .See Table XXIV of the Addendum for updated results.

TABLE XVII .
Results for operators with two derivatives based on the RI ′ -SMOM scheme, obtained by means of the fixed-scale method without the perturbative subtraction of lattice artifacts.The chiral extrapolation has been performed globally.The first number in parentheses gives the statistical error, while the second number is an estimate of the systematic uncertainty.All values refer to the MS scheme at the scale µ 2 0 = 4 GeV 2 .
[86]larly, the five-loop MS anomalous dimension of the quark field in the Landau gauge can be found from Ref.[86].In our paper we used only the four-loop result.The five-loop coefficent reads