Non-perturbative renormalization of three-quark operators

High luminosity accelerators have greatly increased the interest in semi-exclusive and exclusive reactions involving nucleons. The relevant theoretical information is contained in the nucleon wavefunction and can be parametrized by moments of the nucleon distribution amplitudes, which in turn are linked to matrix elements of local three-quark operators. These can be calculated from first principles in lattice QCD. Defining an RI-MOM renormalization scheme, we renormalize three-quark operators corresponding to low moments non-perturbatively and take special care of the operator mixing. After performing a scheme matching and a conversion of the renormalization scale we quote our final results in the MSbar scheme at mu=2 GeV.


Introduction
Distribution amplitudes play an essential role in the investigation of the internal nuclear structure. Exclusive high-energy processes can be factorized into hard and soft subprocesses, where the hard subprocess can be evaluated perturbatively and is characteristic for the reaction in question. The soft subprocess, described by the nucleon distribution amplitude Φ(x i ), contains the information about the distribution of the three valence quark momentum fractions x i inside the nucleon [1,2,3,4,5].
The great interest in this quantity stems from its importance for, e.g., the calculation of the electromagnetic form factors of the nucleon and their scaling behavior. These form factors describe a nucleon absorbing a virtual photon of squared momentum −Q 2 while remaining intact. According to [2], for Q 2 → ∞, the magnetic form factor of the nucleon can be written as a convolution of three amplitudes: first, the distribution amplitude Φ for finding the nucleon in the valence state with the three quarks having definite momentum fractions x i , second, the hard scattering kernel T H , which describes one of the three quarks absorbing the photon, and finally the complex conjugate of Φ that gives the amplitude for the outgoing quarks to form a nucleon again: Here [dx] ≡ dx 1 dx 2 dx 3 δ(1 − i x i ),Q x ≡ min i (x i Q) and m N denotes the nucleon mass.
Distribution amplitudes are genuinely non-perturbative quantities. Hence they are inaccessible for perturbation theory and must be calculated by other means, e.g., by lattice QCD. After a pioneering study in the late 1980s [6] only recently the first quantitative results from lattice QCD have been published [7,8,9]. After performing an expansion near the lightcone, moments of the nucleon distribution amplitudes are expressed in terms of matrix elements of local three-quark operators that are evaluated between a nucleon state and the vacuum. These three-quark operators typically consist of a combination of covariant derivatives acting on the three quark fields, which are located at a common space-time coordinate x. Furthermore, the operators for the nucleon are color singlets and have isospin 1/2.
As the three-quark operators pick up radiative corrections and are subject to mixing with other operators, their renormalization is a vital ingredient for any lattice calculation. In a previous paper [10] we have derived irreducible multiplets of three-quark operators with respect to the spinorial hypercubic group H(4) and have discussed the mixing properties of these operators in detail. The present work will focus on the nonperturbative renormalization of these three-quark operators in the isospin-1/2 sector.
The paper is organized as follows. In the first section, we will introduce an RI-MOM renormalization scheme which will be implemented in our lattice calculations. Then we will explain how to perform a perturbative scheme matching to MS and derive the anomalous dimensions of the operators in question. The following two sections will focus on a discussion of our results for the renormalization matrices, which are finally quoted at µ = 2 GeV in the MS scheme. In the final section we will demonstrate how to renormalize moments of the nucleon distribution amplitude and present the results of consistency checks between the renormalized moments.

Lattice Renormalization
Our main aim is the derivation of non-perturbative renormalization coefficients for three-quark operators in lattice QCD and their application to the renormalization of moments of the nucleon distribution amplitude. As discussed in the introduction, a subsequent convolution with the hard scattering kernel leads then, e.g., to estimates for the electromagnetic form factors of the nucleon. Therefore the renormalization of the hard scattering kernel and the three-quark operators inside the nucleon distribution amplitudes must be carried out in a consistent way, i.e., the same renormalization scheme has to be applied. As the perturbative results for the hard subprocess are usually given in the MS scheme, the distribution amplitudes must also be renormalized in this scheme.
On the lattice, however, it is not possible to implement the MS scheme. Therefore one introduces the RI-MOM renormalization scheme [11]. It is applicable both on the lattice and in the continuum. Then one first renormalizes the matrix elements on the lattice within this scheme. Afterwards one applies continuum perturbation theory to calculate a matching function between both schemes and extracts non-perturbatively renormalized lattice operators in the MS scheme. In this section we will discuss the setup and implementation of the lattice renormalization.

The Three-Quark Operators
We are interested in three-quark operators for the nucleon. Choosing the flavors u, u and d for definiteness, their general form is with analogous definitions for O udu i (x) and O duu i (x). After projection onto isospin 1/2 one ends up with operators for the proton, and subsequent exchange of the u and d fields leads to neutron operators. The color indices c ′ j will be suppressed when not required. The coefficient tensor T (i) determines the symmetry structure of the spinor indices α ′ , β ′ , γ ′ and of the space-time indices µ j , ν j , λ j of the covariant derivatives D contributing to the operator O i .
In [10] we have derived multiplets of three-quark operators that transform irreducibly under the spinorial symmetry group of the hypercubic lattice H(4). These operators reduce the problem of mixing under renormalization to mixing among multiplets belonging to equivalent representations of H(4). Hence they define our choice for the coefficient tensors T (i) . Appendix B of [10] contains all linearly independent sets of potentially mixing isospin-1/2 three-quark operators with leading twist and up to two derivatives. They provide the basis for the (non-perturbative) renormalization. Operators with total derivatives ∂ µ are automatically taken into account due to the identity ∂ µ (f gh) = (D µ f )gh + f (D µ g)h + f g(D µ h) for any color-singlet operator made of the three quark fields f , g and h. Note however that this continuum relation holds only up to discretization errors on the lattice.
Let us finally give an example of a typical three-quark operator belonging to an irreducibly transforming multiplet. We have, e.g., the following operator with two derivatives: The superscript MA (standing for "mixed antisymmetric") indicates that the operator has isospin 1/2 while the subscript f g means that the derivatives act on the first and second quark, compare [10]. Pauli matrices σ µ are used to contract the covariant derivatives. The curly braces denote independent total symmetrization in the dotted and undotted indices of the Weyl representation.

Calculational Method
We introduce a correlation function for the non-perturbative renormalization of the three-quark operators by contracting the operators with three external quark sources, namely u(z 1 ), u(z 2 ) and d(z 3 ): We have two u-quark lines and one d-line running from the three space-time coordinates z i into a vertex at x. There they are connected according to the coefficient tensor We proceed in analogy to the case of quark-antiquark operators [11,12]. We impose fixed momentum on the external quark lines and evaluate the correlation function in momentum space. The result is a four-point function G(p 1 , p 2 , p 3 ) αβγ that depends on the momentum of the three external quark lines and carries three spinor indices: This four-point function is obtained as the ensemble average of quark field contractions on the individual gauge configurations. Let us denote such a contraction on a single configuration by the brackets [.
The most expensive step in the calculation is the evaluation and symmetrization of the spin-color combinations which is, naively speaking, due to the existence of an additional quark and antiquark field with 12 spin-color indices, by a factor of 12×12 = 144 more expensive than for mesonic operators. So, special care has to be taken when implementing these contractions. As the correlation function is not gauge invariant, all configurations are gauge fixed to Landau gauge.
In a final step we amputate the external quark lines of the four-point function G (i) to arrive at the three-quark vertex Γ (i) : where the quark propagators are defined by Note that we can reuse the quantity K(x, p) introduced in (7) to calculate the propagator S(p): In the following section we will define our renormalization scheme based on the vertex Γ (i) .

An RI-MOM Renormalization Scheme
We want to introduce a renormalization scheme that is applicable on the lattice and in the continuum. For quark-antiquark operators such a scheme has been proposed in [11] and is widely known as the RI-MOM scheme. In [13] it was also used for the renormalization of proton decay matrix elements. To study the mixing of our threequark operators with up to two derivatives, we have slightly modified this approach. In the following, we will set up our modified RI-MOM renormalization scheme. As, in general, mixing is a central issue of the renormalization, we will consider mixing matrices explicitly from the very beginning.
The renormalized counterpart of a general regularized operator O i is given by Here Z ij is the renormalization matrix and operator mixing shows up in non-vanishing off-diagonal elements of Z.
We will define the renormalization matrix Z mRI for three-quark operators by projections of the lattice-regularized three-quark vertex Γ introduced in eq. (10). In the following we will distinguish between the tree-level, the lattice regularized (lattice spacing a) and the renormalized vertices Γ tree i , Γ latt i and Γ mRI i , respectively. Let us now introduce a set of projectors P k in spinor space that fulfill the following orthogonality condition with the tree-level vertices: At some renormalization scale µ fixed by the mean squares of the three external quark momenta we then require the renormalized three-quark vertex to fulfill the same equation. This yields the renormalization condition With Γ mRI i = Z Γ,mRI ij Γ latt j we can introduce the auxiliary variable Z Γ,mRI ij for the renormalization of the vertex: In any scheme Z Γ is related to the renormalization matrix of the three-quark operators O by the quark field renormalization Z q . To compensate for the amputated quark legs one needs a factor of Z 1/2 q for each of them: This fixes the three-quark operator renormalization matrix in the RI-MOM scheme: As usual we determine the factors Z q in the RI ′ scheme [11]. The peculiarity of the mRI scheme lies in the special definition of the projectors for the three-quark vertex Γ and will be explained in the next subsection.

Choice of the Projectors
In order to determine the RI-MOM renormalization matrix Z we have introduced but not yet defined a set of projectors P j . The only restriction is given by eq. (14).
We now turn to vector notation. Γ i is a tensor of rank three in spinor-space and can be interpreted as a vector v Γ i of dimension 4 3 . Then we can interpret the projectors P k as orthogonal projections onto vectors v P k of the same dimension: where The task is now to construct a set of vectors v P k that fulfills the normalization condition eq. (14), which reads in vector notation We choose the vectors v P k as follows. We start with an auxiliary vector and project it onto the orthogonal complement of the space spanned by the vectors v Γ tree j , j = k. This results in an altered vector v ′′ P k . Taking care of the normalization we finally define By now all constituents of the mRI renormalization scheme are defined. We summarize the method by rewriting eq. (18) in vector notation: Note that the above renormalization condition will in general depend on the geometry of the external momenta, i.e., the angles between the four-momenta p k . In the end this dependence will be cancelled by the scheme matching. We want to stress furthermore that, due to its general structure, the method is not limited to the case of three-quark operators and four-point functions discussed in this paper, but is applicable to the general class of n-point functions.

Perturbative Calculations
In the introduction we have emphasized the importance of renormalized three-quark operators for nucleon distribution amplitudes. To calculate observables, both the distribution amplitude and the hard scattering kernel must be given in the same renormalization scheme. In the following we will explain, how a matching of the RI-MOM scheme to MS can be achieved with the help of continuum perturbation theory. Moreover we will study the dependence of the renormalization coefficients on the renormalization scale µ. As for the lattice computation, also here all calculations have to be carried out in the Landau gauge.

Scheme Matching to MS
Let us start by writing eq. (13) explicitly in both renormalization schemes: with O mRI i and O MS i denoting the renormalized operators in the modified RI and the MS scheme, respectively. If we introduce the scheme matching matrix we get a relation between the operators renormalized in the MS and mRI schemes: In the following we derive the matching functions with dimensional regularization in continuum perturbation theory.
To this end we proceed with a one-loop perturbative expansion of Z MS and Z mRI . Up to O(ǫ) the renormalization matrix reads in the mRI scheme: Here we have used the renormalized strong coupling α s (µ) = g R (µ) 2 /4π and adopted the following conventions for the dimensional regularization: In the numerical evaluation of α s we use Λ MS = 261 MeV [14]. Comparing eq. (28) with the analogous expression for Z MS and noting that the conversion between both must be finite, the scheme matching matrix in first order becomes This means that Z MS←mRI ij can be derived from the perturbative expansion of Z mRI ij alone. In the following section we will discuss this in more detail.

Determination of Z MS←mRI
According to eq. (17) the renormalization of the three-quark operators consists of four parts: The scheme matching can be performed independently for each renormalization factor so that In two-loop order the matching of the quark field renormalization reads [15]: where ξ is the covariant gauge parameter, compare Appendix A. The scheme matching for the renormalization coefficient Z Γ of the three-quark vertex will be determined in one-loop order. In analogy to eq. (31) we have To evaluate this expression, the renormalization matrix for the three-quark vertex Γ is needed in the mRI scheme. Therefore we perform a perturbative expansion of the dimensionally regularized three-quark vertices: If we apply the projectors introduced in eq. (14) and make use of their linearity, we find Comparing with eq. (28) reveals the identities which have to be evaluated at µ 2 = (p 2 1 + p 2 2 + p 2 3 )/3 for the momentum geometries used in the simulations. Inserting this into eq. (35) yields the result for the scheme matching matrix of the vertex in first order: Together with eqs. (34) and (33) this determines our scheme matching for the threequark operator.
The perturbative calculation of Γ dim i,0 (µ, p k ) for the different three-quark operators is carried out in Euclidean space-time with off-shell quarks and gluons in Landau gauge ξ = 0. It results in lengthy expressions and not all occurring integrals over the Feynman parameters can be solved analytically in closed form. The final evaluation of the integrals over the Feynman parameters and the construction of the projectors P j as well as the evaluation of v P j , v Γ dim i,0 (µ,p k ) were performed with Mathematica for all required momentum combinations.
Note that it is important to exercise care when evaluating the terms of order O(ǫ 0 ) in the three-quark vertex that define Γ dim i,0 . Generally speaking one has to keep track of all possible terms in the Dirac structures proportional to ǫ that multiply a 1 ǫ divergence, since these produce additional contributions to the scheme matching matrix. We have used the following strategy to treat the continuation to d dimensions. In a first step we have written our irreducible three-quark operators in four dimensions as linear combinations of the following basis operators: As usual α, β, γ and τ denote spinor indices, λ i , µ i , ν i as well as ρ and µ are space-time indices. Note that the operator W is equal to the operator V up to the position of the down quark. For the special case of operators without derivatives we have also used to access the operators of sub-leading twist.
Then we have rewritten the three-quark vertices belonging to the above operator basis in terms of dimensionally regularized loop integrals. The corresponding Feynman diagrams consist of three quark lines and one gluon exchange, leading to three strings of gamma matrices in the associated amplitudes. In the vertex two of these strings get contracted due to the presence of the (C . . . ) αβ structure. We can evaluate the remaining contractions of space-time indices using the d-dimensional Dirac algebra with an anticommuting γ 5 and the relation This allows us to identify all contributions that are proportional to ǫ 0 . Once these are determined we construct the regularized vertices of the irreducible three-quark operators from the linear combinations of the A, V , W and U operators that were derived at the very beginning. Finally we evaluate the projections of eq. (39) resulting in the desired scheme matching matrices Z Γ,MS←mRI .

Renormalization Group Behavior
Knowing the scaling behavior of Z provides a valuable consistency check of the results. To this end one can compare lattice results that were derived at different renormalization scales with the perturbatively expected scaling. This will be done in subsection 4.2.
Generally, the scaling behavior of a quantity is described by the renormalization group equation (RGE) and the related beta-and gamma-functions Both β and the anomalous dimension γ can be written as an expansion in the strong coupling: In the MS scheme, α s as well as the beta function are known to high order in perturbative QCD. The gamma function is operator dependent and in most cases not known to the same accuracy.
Let us focus on the behavior of the renormalization matrix for the three-quark operators under a change of the renormalization scale. It follows from the RGE that a renormalization matrix can be converted from one scale µ to any other scaleμ with the help of a scaling matrix ∆Z MS : The matrix ∆Z MS (µ) again depends on the beta-function and the operator-specific gamma-function. We derive the scaling matrix for the three-quark operators by independently converting the renormalization matrices of the three-quark vertex and the three quark fields, cf. eq. (17). As the behavior of the quark fields is known to higher accuracy, we hope to get a better description for the three-quark operator renormalization by treating this contribution, just as the quark field scheme matching, to order α 2 s .
Hence we apply the expression (A.4) in Appendix A.2 to convert the wave-function renormalization Z q together with its anomalous dimension from Appendix A.3. A leading order formula is used for the scaling of the three-quark vertex renormalization matrix Z Γ,MS (µ): We will use this to rescale our results obtained at a set of different µs to the final scalẽ µ = 2 GeV.
We still need the gamma function γ Γ of the three-quark vertex. Omitting contributions that vanish for ǫ → 0, eq. (36) can be cast in the form where X 2 is some momentum square, µ the scale of the α s expansion and C j is a finite term. We thus find for the renormalization matrix in the MS scheme Using eq. (44) and g R (µ) 2 = g 2 µ −ǫ +O(g 4 R ) yields the anomalous dimension matrices: The results for the anomalous dimensions are summarized in Appendix B.
The whole procedure can only be expected to work in a "window", where the renormalization scale is large enough for perturbation theory to be a good approximation and small enough so that the lattice cutoff is sufficiently far away. Hence we expect the condition to restrict µ to a reasonable range.

Details of the Lattice Setup
We work with non-perturbatively O(a) improved clover-Wilson fermions and the plaquette gauge action. We use gauge configurations generated by the QCDSF/UKQCD collaborations with two dynamical flavors of sea quarks, n f = 2. The gauge is fixed to Landau gauge by minimizing the functional with Table 1 The used gauge configurations of the QCDSF/UKQCD collaborations. The first column shows the available lattice couplings β, the second column the lattice volume and the following columns summarizes the hopping parameters κ.  Table 2 Critical hopping parameters κ c and inverse lattice spacings in units of the Sommer parameter r 0 extrapolated to the chiral limit. For any gauge link U µ (x) connecting the sites x and x +μ one can construct gauge equivalent links U G µ (x) by applying SU (3) transformations G(x). Landau gauge is realized by the gauge configuration U G µ (x) that minimizes the functional F G [U] [16,17].
We have performed computations on the ensembles of gauge field configurations listed in Table 1. Table 2 provides further information on the chiral limit of these lattices. Our approach uses momentum sources [12] which improves the statistics significantly, as it "averages" over all space-time coordinates and thus is superior to point-source methods. Due to this benefit an order of ten gauge configurations per ensemble provide sufficiently high statistics for our purposes. Anyhow, the statistical error will turn out to be negligible compared to the systematic error related to the perturbative scheme matching and scale conversion.
All calculations for the renormalization coefficients were performed on the Regensburg QCDOC machine with partitions of up to 128 nodes. The code implementation was done in C++ using QDP++, the SciDAC data-parallel programming interface. We have taken over the lattice action and an even-odd-preconditioned stabilized bi-conjugate gradient method for the inversion of the Dirac-operator from the Chroma software library [18,19].
For the discretized versions of the first and second covariant derivatives acting on a fermion field inside a three-quark operator we have used To minimize discretization effects, we choose the three external four-momenta p i of the quarks close to the diagonals of the Brillouin zone. Varying the squares of the momenta allows us to derive the renormalization coefficients at several different values of the renormalization scale µ. The Sommer parameter r 0 = 0.467 fm [14] is used to set the scale.
At the end we have to perform a chiral extrapolation of the results. Therefore, we compute the amputated three-quark vertices Γ for fixed lattice size and fixed coupling β at typically three different values of the hopping-parameter κ. For each κ value we then determine the "renormalization matrix"Z mRI ij (µ, κ). The final mRI renormalization coefficients Z mRI ij (µ) are derived from a linear extrapolation in 1/κ to the critical value 1/κ c . A typical chiral extrapolation is depicted in Fig. 1 and shows a reasonably linear behavior. Fig. 1. A typical chiral extrapolation for a renormalization coefficient of an operator with two covariant derivatives. The curve shows a linear fit in 1/κ at a fixed scale µ in the RI-MOM scheme, the cross marks the chiral limit.

Data Analysis
The Jackknife method is used to estimate the statistical errors for all elements of the matricesZ mRI ij (µ, κ). From the uncertainty in the χ 2 fit to the chiral limit we then determine the statistical error for each coefficient of the final result Z mRI ij (µ).
Due to its definition, Z mRI ij (µ) may also contain imaginary parts before the scheme matching to MS. This can be seen best from eq. (49). Applying the projector P j to the right-hand-side of this equation yields the inverse of the mRI renormalization matrix. It is obvious that, apart from the real terms resulting from the expressions proportional to Γ tree i , we will also end up with a term proportional to P j C i from the finite contribution. Written in vector notation this projection reads v P j , v C i and this scalar product is not necessarily real. So finite contributions like C i may render the mRI renormalization matrix complex. The scheme matching to MS should cancel these terms. However, since we perform a one-loop scheme matching only, we cannot expect a complete cancellation of the imaginary parts. In practice it turns out that the residual imaginary parts are reasonably small and may be used as an independent estimate of the systematic uncertainty induced by the scheme matching. Hence we will quote the real part of the renormalization matrices Z MS ij (µ) as our results.
The systematic error is estimated via the scaling behavior of the renormalization matrices. To this end we perform a lattice renormalization at ten different renormalization scales µ. We exemplify our procedure with Fig. 2, where we have plotted the scaling of one diagonal and one off-diagonal coefficient of the renormalization matrix for threequark operators with two derivatives in the irreducible representation τ 12 2 . The general behavior is typical for all representations: While the diagonal coefficient is of order In the investigated range µ 2 = 3 GeV 2 − 100 GeV 2 the renormalization coefficients vary by roughly 40%. The scheme matching to MS slightly shifts both mRI curves to lower absolute values (dash-dotted curve). In the third step we convert all results to the same renormalization scaleμ = 2 GeV (solid curve). This leads to a manifest flattening of the graph. Ideally the curve would be a constant now, but as already stated neither the lattice, nor the one-loop perturbative approach is free of systematic uncertainties in the whole µ 2 range. At low scales we expect the perturbative expansion to break down, as the coupling becomes too strong for our one-loop approach to hold. This behavior can be observed in the region µ 2 < 10 GeV 2 . On the other hand we also have to expect cutoff effects from the lattice calculation at large scales. Also this feature can be found in the figure. In the regime of µ 2 = 10 GeV 2 − 40 GeV 2 the results appear almost flat, which may be interpreted as a scaling window. Note that although operators with a different number of derivatives may come along with scaling windows shifted to slightly higher or lower values of µ 2 , all scaling windows seem to overlap in the quoted region. This indicates that the chosen approach for the non-perturbative renormalization and perturbative scheme matching is well justified.
Each renormalization coefficient converted to the scaleμ = 2 GeV can be reasonably described by the three-parameter fitting formula where µ stands for the original renormalization scale before the conversion. This formula incorporates O(a 2 ) effects as well as potential further logarithms in µ 2 stemming from higher order corrections in the scheme matching and scale conversion. Our final result is read off at µ 2 = 20 GeV 2 . As our estimate of the systematic error we take the maximum of the differences with the values at 10 GeV 2 and 40 GeV 2 .
For the two renormalization coefficients plotted in Fig. 2 one finds in this way where the first error is the statistical and the second the systematic uncertainty. As expected for our one-loop calculation, the systematic uncertainty is larger than the statistical error and exceeds it by a typical factor of O(10) to O(100). It would be interesting to test our estimates for the systematic error by comparing them with future higher order calculations of the perturbative expressions.
We have also investigated the dependence of the renormalization matrices on the geometry of the external quark momenta p i . To this end we have calculated the Z matrices for modified angles between the three quark lines running into the vertex. In Fig. 3 we show a typical result for the same renormalization coefficients as in Fig. 2. We find that up to roughly 30 GeV 2 even in the RI-MOM scheme no significant difference between both geometries occurs, whereas at still larger scales small deviations are observed.
in good agreement with the values in eq. (58).

Results for Z MS (2 GeV)
Let us now present the renormalization matrices in the MS scheme at µ = 2 GeV for the different H(4) irreducible representations τ 4 1 , τ 4 2 , τ 8 , τ 12 1 , τ 12 2 of isospin-1/2 threequark operators. For these representations we use the same notation as in Ref. [10]. The superscript denotes the dimension of the representation, with the underscore indicating its spinorial nature. The subscript distinguishes inequivalent representations of the same dimension.
We will denote the statistical error by E st and the systematic error by E sy . As mentioned, the latter is estimated by comparing the Z values at 20 GeV 2 with those at 10 GeV 2 and 40 GeV 2 . Thus the final result reads As the statistical errors are almost two orders of magnitude smaller than the systematic errors, we will quote the statistical errors only in two explanatory examples and drop them for the rest of the paper.
Another source of uncertainty results from the error of Λ MS . We use Λ MS = 261 MeV from Ref. [14]. Adding the two given errors we obtain a combined error of 43 MeV. While the resulting uncertainty is considerably smaller than E sy for operators with zero and one derivative it becomes comparable with E sy in the case of two derivatives.
The mixing multiplets can be read off from Table 3: Multiplets of the same representation and same dimension can mix under renormalization. Lower-dimensional operators can also mix into higher dimensional ones of the same representation by powers of the inverse lattice spacing 1 a . We will summarize these 1 a and 1 a 2 admixtures -wherever present -in the last columns of the related renormalization matrix (cf. the operator basis in Appendix C).
We find that the results are similar for different values of the lattice coupling β and essentially identical at different lattice sizes, which again demonstrates the consistency of the approach. For reasons of better readability, we will only discuss the result for our finest lattice β = 5.40, L 3 × T = 24 3 × 48 in this section, while more details are given in Appendix C. There we also give the operator bases explicitly.

Representation τ 4 1 (Zero Derivatives)
The following renormalization matrix belongs to the three-quark operators without derivatives in the irreducible representation τ 4 1 . There are two multiplets that mix with Table 3 Irreducibly transforming multiplets of three-quark operators with isospin 1/2 sorted by their mass dimension. The subscripts f , g and h indicate that the covariant derivative(s) act on the first, second or third quark, respectively, cf. [10]. dimension 9/2 dimension 11/2 dimension 13/2 The diagonal elements of Z are smaller than one and for both operators of almost equal size. The mixing off-diagonal elements are both negative and amount to roughly four and one per cent of the diagonal coefficients. Furthermore we see that the statistical error is of relative order 10 −4 which renders it negligible compared to the systematic error. The latter one is two orders of magnitude larger and hence will be the dominating source of uncertainty in any renormalization of matrix elements on the lattice.

1 (Zero Derivatives)
There is only one multiplet of three-quark operators without derivatives belonging to the irreducible representation τ 12 1 . Therefore the renormalization matrix becomes onedimensional. Our result reads: Again, the diagonal element is smaller than one and the statistical error is of the order 10 −4 so that the systematic error dominates.

Representation τ 8 (One Derivative)
We carry on with three-quark operators with one covariant derivative. Since the statistical and systematic errors are of similar order of magnitude as for the above operators without derivatives, we do not quote them here anymore for reasons of better readability, but refer to the summary in Appendix C. We obtain Z = 1.1260 .
The diagonal component of the renormalization matrix for operators with one covariant derivative is larger than those for operators without derivatives.

1 (One Derivative)
For leading-twist operators with one covariant derivative there are three multiplets of the representation τ 12 1 . One lower-dimensional multiplet that belongs to the same irreducible representation and is formed by operators without covariant derivatives can mix with them via one power of the inverse lattice-spacing 1 a , cf.  The coefficients in the last column, which describe the mixing with lower-dimensional operators, are rather small.

2 (One Derivative)
This representation has four mixing multiplets. As for τ 8 there is no mixing with lowerdimensional operators. Therefore these two representations are especially well suited for the evaluation of matrix elements involving three-quark operators with one derivative. We find: The mixing of the second multiplet with the first one is in the realm of ten per cent, while the mixing of the fourth multiplet can almost be neglected. This might be related to the fact that it contains three-quark operators with different chiralities of the quark fields than the other three multiplets.

Representation τ 4 1 (Two Derivatives)
We now proceed with leading-twist operators with two covariant derivatives. These operators belong to five inequivalent irreducible representations of the spinorial hypercubic group H(4). We start with the six multiplets of operators with two derivatives in τ 4 1 , which can mix with the two lower-dimensional multiplets, which contain threequark operators without derivatives. The 1 a 2 -admixture of the latter operators is seen in the last two columns of the renormalization matrix Z: Here the statistical errors are by almost a factor of ten larger than for operators without derivatives. However, the relative systematic errors do not seem to change considerably. We furthermore observe a clear hierarchy with the diagonal elements further increasing relative to three-quark operators with one and without derivatives.

Representation τ 4 2 (Two Derivatives)
This irreducible representation consists of six leading-twist multiplets. As by grouptheoretical arguments mixing with lower-dimensional three-quark operators can be excluded, τ 4 2 is the representation of choice wherever possible when working with threequark operators with two derivatives. The renormalization matrix at 2 GeV reads in the MS scheme:

Representation τ 8 (Two Derivatives)
Again, six multiplets of operators with two derivatives contribute to this irreducible representation. They can mix with the lower-dimensional operators with one derivative of the representation τ 8 : The mixing coefficients in the last column indicate that the 1 a -contribution of the lowerdimensional multiplet might be important in practical applications.

1 (Two Derivatives)
Eight operator multiplets with two derivatives mix with the three multiplets with one derivative and the operator multiplet without derivatives of τ 12 1 . For reasons of better readability we split off the last four columns of the renormalization matrix, which describe the mixing with these lower-dimensional operators, and summarize the coefficients in a separate matrix Z ′ . Then the 1 a mixing with the operators with one derivatives can be read off from the first three columns of Z ′ , the 1 a 2 mixing from its last column: Compared to the diagonal elements around 1.3, the sixth operator mixes into the first with 0.2 and the 1 a 2 mixing contributes with a factor of up to 0.02. When renormalizing matrix elements with these coefficients, the mixing of the lower-dimensional lattice operators can only be neglected, if the matrix elements of these operators are considerably smaller than those containing operators with two derivatives. Whether this is the case or not must be decided in each individual case. As already stated earlier, whenever possible one should anyway make use of the operators of representation τ 4 2 , because these do not suffer from mixing with lower-dimensional operators.

2 (Two Derivatives)
This last representation also contains eight three-quark operators with two derivatives that can mix with the four lower-dimensional operators with one derivative of the same representation. For reasons of better readability we have again split off the last four columns of the renormalization matrix and display the related coefficients in the sepa-rate matrix Z ′ : Also in this case the renormalization matrix has mixing coefficients of a few per cent for the lower-dimensional operators. Whether or not they have to be taken into account depends again on the magnitude of the corresponding matrix elements.

Renormalization of Moments of the Nucleon Distribution Amplitude
Having presented the renormalization matrices in the previous section, the next step is their application to the moments of the nucleon distribution amplitude φ: For a detailed discussion of the calculation of these quantities in lattice QCD, the notation used and the physical interpretation of the results we refer to [8,9]. Here we focus on the behavior under renormalization and quote numbers only to exemplify the involved orders of magnitude.

Zeroth Moment
The zeroth moment is linked to a matrix element of three-quark operators without derivatives that is proportional to f N , the normalization constant of the nucleon wave function: Here N denotes the nucleon spinor. The definition of O 000 A,0 can be found in [9]. For the renormalization of the matrix element only the multiplet and the irreducible representation it belongs to is relevant. The first spinor component reads in our MA-isospin operator-basis: This operator belongs to the irreducible representation τ (67) Table 4.

Next-to-Leading Twist Constants λ 1 , λ 2
For three-quark operators without derivatives we have also computed the renormalization matrix for next-to-leading twist. This enables us to renormalize the constants λ 1 and λ 2 , which describe the coupling of the nucleon to two different interpolating fields used in QCD sum rules. Since we work in Euclidean space we have: We perform the following discussion in detail to clarify the treatment of mixing matrix elements and the required change of basis. In a first step we relate the operators sandwiched between the nucleon and vacuum to our MA-isospin basis: Hence λ 1 renormalizes like −8 O we have Substituting back the bare λs yields finally the desired relation between the bare and renormalized values:

First Moments
The first moments of the nucleon distribution amplitude can be derived from matrix elements of three-quark operators with one covariant derivative: The operators O lmn A,1 , where l + m + n = 1, are related to the MA operators of the representation τ , The matrix elements of these operators yield the bare values for f N φ 100 , f N φ 010 and f N φ 001 , as eq. (73) shows. Using the definitions  (20) we obtain the following relation to the renormalized first moments of the nucleon distribution amplitude: Upon using f ren N = Z(τ 12 1 ) f N we get: Due to eqs. (63) and (64) the moments must comply with While the bare values do not fulfill this equation (the sum equals 0.84), the sum of the renormalized values is found to be 0.97. Within errors this is in agreement with the constraint.

Second Moments
Analogously we can renormalize the second moments of the nucleon distribution amplitude. A typical matrix element is given by Again we investigate the multiplets and the representation to which the operators O lmn 2 , l + m + n = 2, belong by rewriting them in terms of our isospin mixed-antisymmetric basis. For the fourth spinor component we have, e.g., With the new observables we find the renormalized quantities Here, the renormalized values must fulfill four constraints: Comparing the results in Table 4 reveals good agreement within errors also for these quantities. This is an encouraging result, which demonstrates the consistency of the applied renormalization. It also makes us confident that we did not severely underestimate our errors.

Summary and Conclusions
In this paper we have set up a lattice renormalization scheme for three-quark operators based on the RI-MOM approach. The introduction of a suitable set of projectors fa-cilitates the definition and practical evaluation of a renormalization matrix in our mRI scheme. In a second step we have performed a scheme matching to MS and a scale conversion to 2 GeV based on one-loop continuum perturbation theory. It was found that the statistical errors are negligible compared to the systematic uncertainties. Finally we have explained how to apply our results to the renormalization of low moments of the nucleon distribution amplitude.

A.1 The Action
In this appendix we summarize the conventions that we have adopted for the perturbative evaluation of the three-quark vertex Γ in one-loop order. We use the Euclidean action Here ψ denotes the fermion fields, u the ghosts, A a µ the gluon field and λ a stands for the Gell-Mann matrices fulfilling [λ a , λ b ] = if abc λ c . The propagators in momentum space are given by: We work in the chiral limit m = 0 and use the Landau gauge ξ = 0.

A.2 The Scaling Function
A three-loop formula for the function ∆Z M S (µ) is given by: with the abbreviationsβ i = β i /β 0 andγ i = γ i /(2β 0 ).

A.3 The Quark Field Anomalous Dimension
Finally we give the anomalous dimension of the quark field so that we can extract the scaling behavior from eq. (A.4). It reads [20]: ). Here

B Anomalous Dimensions
In the following we summarize the gauge invariant anomalous dimensions γ for the three-quark operators. In the one-loop perturbative expansion they are identical to the anomalous dimension γ Γ , eq. (51), of the three-quark vertex in Landau gauge: The operator basis is taken over from the linearly independent, isospin-1/2 H(4) irreducibly transforming multiplets of three-quark operators introduced in Appendix B of [10]. We have checked that the eigenvalues of our anomalous dimension matrices agree with those presented by Peskin in [21].

B.1 Operators without Derivatives
We start with three-quark operators without covariant derivatives. For the irreducible representation τ 4 1 we choose the operator basis MA 3 .
To first order we find the anomalous dimension matrix The operator basis for the representation τ . Then

B.2 Operators with One Derivative
In the following we summarize the anomalous dimensions for three-quark operators with one covariant derivative. The representation τ 8 has the operator basis form a basis and have the anomalous dimension matrix has the one-loop anomalous dimension The same anomalous dimension can be shown to apply to the second twelve-dimensional representation τ DD12 2 and its operator basis

C The Renormalization Matrices for Three-Quark Operators
In this appendix we present the renormalization matrix Z(2 GeV) and the error matrices E sy . As the statistical errors are much smaller than the systematic uncertainties we do not quote them here. Note that the error due to the error of Λ MS leads to an additional non-negligible uncertainty (similar in size to E sy ) for operators with two derivatives. We also list the operator bases for all irreducible representations. For the notation of the three-quark operators compare again [10]. The renormalized operators are related to the bare lattice operators by We have results not only for the two lattices presented in the following, but also for all other lattices in Table 1. However, in order to keep the paper at a reasonable length we restrict ourselves to the largest lattices (24 3 × 48) at β = 5.29 and β = 5.40. 4 1 This irreducible representation contains two mixing multiplets. The renormalization matrix Z ij is given in the following operator basis: MA 3 .

(C.2)
We now present our chirally extrapolated results. Only one operator multiplet belongs to τ 12 1 . Hence, there is no mixing present and the operator basis is given by MA 7 .
(C.5) In leading twist there is also no mixing for this representation. We take the basis There are four mixing multiplets and we take We work in the following operator basis: Here, mixing with lower-dimensional operators occurs:

2
This is the only irreducible representation of leading-twist operators with two derivatives that is not subject to mixing with lower-dimensional operators on the lattice: Here, twelve multiplets of operators mix with each other under renormalization. Four of them have lower dimension: MA 7 .
We have split off the last four columns of the renormalization matrix, which describe the mixing with the lower-dimensional operators O 9 ,...,O 12 , and display the related coefficients in a separate matrix Z ′ . β = 5.29, lattice size:    Finally, we have the representation τ