Moments of Ioffe time parton distribution functions from non-local matrix elements

We examine the relation of moments of parton distribution functions to matrix elements of non-local operators computed in lattice quantum chromodynamics. We argue that after the continuum limit is taken, these non-local matrix elements give access to moments that are finite and can be matched to those defined in the $\overline{MS}$ scheme. We demonstrate this fact with a numerical computation of moments through non-local matrix elements in the quenched approximation and we find that these moments are in excellent agreement with the moments obtained from direct computations of local twist-2 matrix elements in the quenched approximation.

It should be noted that the same matrix element is used to define the collinear PDFs by taking z to be a separation along the light cone. In particular, with z = (0, z − , 0, 0) and α = + we obtain with the second term dropping out because z does not have a "+" component and noting that z 2 = 0 on the light cone. Given that M p (zp, z 2 ) is Lorentz invariant, we can compute it with any convenient choice of z and α, however, if z is not on the light cone, the limit z 2 → 0 has to be taken. This limit is non-trivial as there exists a logarithmic singularity at z 2 = 0. For that reason the PDFs are defined through factorization of this short distance singularity which results in the scale dependent PDFs [17][18][19]. Furthermore, it is well known that in the limit of z 2 = 0 only the twist-2 contribution survives. With the above discussion it is easy to see that a non-perturbative computation of the collinear PDFs should start with the computation of the invariant function M p (zp, z 2 ) from which one can obtain the twist-2 contribution in the limit z 2 → 0. Before continuing further, it should be noted that the Lorentz invariant quantity zp is called Ioffe time ν in the literature [20,21]. A convenient choice is z = (0, 0, 0, z 3 ), α in the temporal direction i.e. α = 0, and the hadron momentum p = (p 0 , 0, 0, p). In this case the z α -part drops out isolating the function M p (ν, z 2 ) we seek to compute [33]. Considering a time local matrix element as above also allows for its computation in Euclidean space using lattice QCD. As it is well known [22], in this case the matrix element computed in Euclidean space is the same as the Minkowski space counterpart. Therefore, one can analyze the properties of the above matrix element using the Minkowski metric. Unfortunately, taking z to be space-like, UV singularities arise from the gauge link self energy and end points [18,19,23,24]. It has been shown, that these UV singularities can be factorized in a multiplicative factor and can be renormalized away. A particularly practical proposal for removing these singularities is by considering the ratio which exactly cancels all these UV singularities leaving behind a finite reduced function for which the regulator can be removed and therefore, when computed in lattice QCD it can be extrapolated to the continuum limit at fixed ν and z 2 . This approach was discussed and tested in [16,25], where it was indeed shown that these UV singularities are absent from the numerical data for M(ν, z 2 ). Furthermore, with the UV singularities canceled in the ratio the only remaining singularity at z 2 = 0 is that associated with collinear divergences and therefore the operator product expansion (OPE) can be used on M(ν, z 2 ) resulting in a factorization into a PDF and a perturbatively computable coefficient function. Therefore we can write where µ is the factorization scale in a particular scheme such as M S, α s (µ 2 ) is the strong coupling constant and Q(ν, µ) is the Ioffe time PDF in that scheme. Furthermore, there are additional polynomial corrections to the leading order expression that vanish in the limit of z 2 = 0. These corrections are not the same for M(ν, z 2 ) and M p (ν, z 2 ). As it was shown in [16,25], the polynomial corrections for the case of the reduced function are small due to cancellation of certain polynomial terms in the numerator and the denominator. In particular, because the reduced function by construction satisfies M(0, z 2 ) = 1, one can see that B k (0) = 0 and 1 0 dα C(α, z 2 µ 2 , α s (µ)) = 1 .
In the subsequent subsections we proceed with further discussion of UV divergences and consider the expansion in to moments of the reduced function M(ν, z 2 ) from where the M S moments of PDFs can be extracted.

A. Comment on UV divergences
The matrix element M p (ν, z 2 ) has UV divergences that need to be renormalized before any discussion of the collinear PDFs begins. However, as we said before, these divergences can be renormalized multiplicatively. Given that there exist many ways to treat these divergences one may want to look at the most general case scenario. Let's assume we adopt a particular regularization scheme (a lattice with lattice spacing a) and a particular renormalization scheme (e.g. RI-MOM [26]). In this case, the renormalized matrix element at a scale µ U V (in the case of multiplicative renormalization) would be related to the bare matrix (denoted by the subscript a) element by The renormalized matrix element can then be extrapolated to the continuum limit as all UV divergences have been removed. It should be noted that the renormalized matrix element depends now on a new scale, the renormalization scale, which takes the place of the lattice spacing. This scale is different from the factorization scale at which the PDFs are defined as one renormalizes UV divergences while the other factorizes collinear divergences present in the PDFs.
In addition, one may define the renormalization group invariant matrix elements (RGI) through the continuum UV invariant scaling functions σ(µ U V ) (see for example [27]) by where σ(µ U V ) is regulator independent, but scheme dependent, and the RGI matrix element is scheme and regulator independent. This RGI matrix elements are the objects that need to be analyzed in order to obtain the collinear PDFs. One may obtain them in a variety of methods (ex. using the RI-MOM scheme [28]), and can discuss their properties without the need for referring to the lattice regulator anymore. Therefore, the ratio is regulator and RG independent and does not have any knowledge of a lattice cut-off that may have been used in order to compute it. Here M RGI p is defined by With this discussion it is clear that the expansion of M(ν, z 2 ) into leading twist parton distribution functions and subleading higher twist contributions that are suppressed by powers of z 2 cannot fail due to lattice artifacts as it can be performed without the use of a lattice regulator.

B. Expansion in moments and power divergences
Although the argument from the last section is conclusive in understanding the issue of power UV divergences in direct computations of PDFs, it is still instructive to see how power divergent mixing present in the lattice regularized moments can be eliminated in the full z-dependent matrix element. For this reason in this section we will discuss the OPE of the reduced function M(ν, z 2 ). In the limit of small z 2 , the matrix elements of the numerator and the denominator can be expanded using OPE in terms of local non-perturbative, renormalized matrix elements and Wilson coefficients. The Wilson coefficients can in principle be computed in perturbation theory. The matrix elements can be computed in any regularization and renormalization scheme we desire. However, if a lattice regulator is used, power divergences may arise. In particular, if a lattice regulator is used then the breaking of rotational symmetry introduces mixing between different spins even at the leading twist. In our analysis, we will first ignore higher twist effects and focus only onto the impact of power divergences due to mixing between operators with different spin at twist 2. The breaking of O(4) symmetry on a Euclidean lattice collapses the infinite number of continuum irreducible representations of the rotation group to a set of finite irreducible representations of the hypercubic group H(4). Operators that belong in a particular H(4) irreducible representation can always be written as a linear combination of the continuum operators of definite spin. As a result, operators of different dimensions can mix under renormalization, resulting in power divergences as the continuum limit is taken. The mixing of a particular operator can be classified by the difference in dimension between it and the operator whose matrix element it is contaminating. Mixing with higher dimensional operators comes with positive powers of the lattice spacing and hence is eliminated in the continuum limit. Therefore, we can ignore this mixing. On the contrary, mixing with operators of the same or lower dimension survives in the continuum (or even worse, diverges). Therefore, special care has to be given to operators of lower or equal dimension.
If we ignore higher twist effects, then the small z 2 expansion of our matrix element would read as where p|O 0α1···α k (k) |p µ are the familiar spin averaged twist-2 matrix elements, and µ is the factorization scale. The above equation is valid if a factorization scheme without power divergences is chosen. Such a scheme is M S. Note that the factor of 1/(2p 0 ) is inherited from the denominator matrix element that defines M(ν, z 2 ), which is also expanded in powers of z. However, because the denominator is a zero momentum matrix element, it does not contain a tower of twist-2 matrix elements which all vanish at zero momentum with the exception of the vector current matrix element which is equal to one in the isovector case that we are considering here. Furthermore, the higher twist effects in the denominator are considered small and are reabsorbed in the O(z 2 ) terms that are omitted. In the previous section we argued that the coefficients B k (ν) in the higher twist expansion are such that B k (0) = 0. It is precisely the reabsorption of the denominator polynomial in the z 2 term that cancels which brings about B k (0) = 0. Therefore, we can explicitly see how part of the higher twist effects is canceled in the M(ν, z 2 ).
We could write down the same equation using the lattice regularized matrix elements as where now a is the lattice spacing andc k are the Wilson coefficients of the lattice expansion. In the above expression, the matrix elements P |O 0α1···α k (k) |P a , which are finite at fixed a, do not have a well defined continuum limit. However, the left hand side in the equation is a well defined quantity in the continuum. Therefore, one should be able to reconcile these two facts by explicitly showing that any power divergent cut-off dependence in the right hand side of the equation cancels. One can show this by considering the relation between the renormalized M S matrix elements with the lattice regularized matrix elements. This reads where B kl are scale independent mixing parameters which subtract the power divergences of the lattice regularized matrix elements, and Z (k) (µ 2 a 2 ) are the M S renormalization factors. Note that the product of c k (z 2 µ 2 )Z (k) (µ 2 a 2 ) does not depend on the scale µ and because the right hand side of Eq. (12) has to be scale independent the Wilson coefficient scale dependence has to cancel the matrix element scale dependence order by order in the expansion. In order to simplify matters, all indices in O 0α1···α k (k) should be considered to be along the same direction (lets say α k = 3 and z a k = z for all k). We can now put the above expression in Eq. (12) and rearrange the sum to make each term proportional to a lattice matrix element. It is clear that each lattice matrix element appears in all M S terms of the same or higher order. It is easy to see that We can now identify the relation of the lattice Wilson coefficients to the M S Wilson coefficients and show that Note that the right hand side is scale independent due to the fact that always the scale independent pair c k (z 2 µ 2 )Z (k) (µ 2 a 2 ) appears. Furthermore one can see that the lattice Wilson coefficients have power divergences. This is not surprising given that we are expanding a UV finite quantity and thus the Wilson coefficients have to absorb the matrix element power divergences. An explicit computation of the lattice Wilson coefficients in perturbation theory will also reveal that fact. At this point it is clear that a more general statement can be made: if the OPE is performed on a UV finite quantity using operators that under the chosen regularization scheme have power divergences, those divergences have to be canceled exactly from corresponding power divergences in the Wilson coefficients. This has been shown in the past in [29]. Therefore, our reduced matrix element M(ν, z 2 ) can be expanded in power divergent lattice operators using OPE provided the corresponding power divergent Wilson coefficients are used. However, it will still be free of power divergences due to the cancellation between divergences in Wilson coefficients and matrix elements.

III. COMPUTATION OF MOMENTS
In this section we discuss the computation of moments of PDFs from M(ν, z 2 ) which can be computed on the lattice and has a well defined continuum limit. Having established that the expansion in moments is well defined in any scheme, we chose to work with M S in this section. First one can further simplify the expression in Eq. (12) by replacing the matrix elements with moments in M S a n (µ). These moments are defined by where [· · · ] sym stands for symmetrization of indices. Inserting this in Eq. (12) we obtain where the product p 3 z 3 has been replaced by the Ioffe time ν. This definition of the moments is derived by the traditional definition a n (µ) = where q(x, µ) is the parton distribution function. Considering the definition of Ioffe time PDFs, we can derive that where a n (µ), is the n-th moment of the parton distribution function. From this expression and Eq. (6) we obtain that if one expands in a Taylor series with respect to ν the reduced function M(ν, z 2 ), the coefficients of this Taylor series expansion are the moments of the PDFs up to a multiplicative constant and up to O(z 2 ) higher twist effects. In other words, Furthermore, Eq. (21) implies that the Wilson coefficients are c n (z 2 µ 2 ) = 1 0 da C(α, z 2 µ 2 , α s (µ))α n .
Since C(α, z 2 µ 2 , α s (µ)) is known analytically [19,24,30] to first order in α s in M S, we can easily compute the Wilson coefficients c n (z 2 µ 2 ) in M S, by simple integration of Eq. (23). The leading order M S expression for C(α, z 2 µ 2 , α s (µ)) is where B(a) is the Altarelli-Parisi kernel and D(α) is given by In the above equations [· · · ] + denotes the "plus prescription". Integrating Eq. (23) one obtains where which are the well known leading order moments of the Altarelli-Parisi kernel, and Here ψ (1) (z) is the polygamma function defined as ψ (1) (z) = d 2 ln Γ(z)/dz 2 with Γ(z) being the Γ-function. With the Wilson coefficients computed we can now obtain the M S moments up to O(α 2 s , z 2 ) directly from the reduced function M(ν, z 2 ) as Note that in order to do so one needs a precise computation of M(ν, z 2 ) in the small ν region at fixed z 2 . This is the region in which lattice computations can easily achieve high precision.  [31]. At low z 2 the perturbative matching seems to work well as indicated by the independence of the moment on z 2 /a 2 .
To illustrate this procedure we take as an example a recent quenched QCD calculation [16] which can be compared with the results from [31] where the moments where obtained through direct computations of the matrix elements of twist-2 operators. In [16], the reduced isovector Ioffe time pseudo-PDF was computed at a fixed coupling β = 6.0 in quenched QCD using the Wilson gauge action and Clover improved valence fermions. The lattice spacing in this computations is 0.093 fm. The same quenched theory was also used in [31] to study the moments of PDFs from direct computations of the corresponding twist-2 nucleon matrix elements, however in this case Wilson fermions were used for the valence quarks. The two calculations have very different systematics and most importantly different discretization errors due to the use of two different valence fermion actions that differ by O(a) effects. Nonetheless, it is instructive to check if the moments computed from the reduced Ioffe time PDF agree within these expected systematic effects with the direct computation. For our comparison the pion mass is set in both cases to m π ≈ 600 MeV.
On the left panel of Fig. 1, we plot the left hand side of Eq. (22). These are the derivatives of M(ν, z 2 ), rescaled by powers of i, at ν = 0. The derivatives of M(ν, z 2 ) are estimated numerically from its real and imaginary parts, using finite difference derivatives with O(ν 4 ) errors, which are of the order of a few percent. The real part contains only even powers in ν while the imaginary only odd. This is taken into account in order to simplify the numerical derivative estimator.
In the right panel of Fig. 1, we plot the lowest two moments of the unpolarized isovector PDF computed using Eq. (30) at scale µ = 4 GeV. These moments are plotted as function of z 2 /a 2 used in their extraction. At this scale the perturbative corrections are of order 10%. At lower scales these corrections become larger as expected and thus the 1-loop matching is expected to break down. Unfortunately, the available data for M(ν, z 2 ) do not have sufficient precision to extract higher moments. As we can see at small z 2 , where the perturbative expansion is expected to work better, the resulting moment is independent of z 2 indicating that the matching formula works well in this region of z 2 . At higher values of z 2 , as expected the perturbative matching breaks down and thus the moment depends on the z 2 used in the computation. Furthermore, the small variation of the derivatives of M(ν, z 2 ) and the extracted moments over a wide range of z 2 indicates that polynomial corrections are indeed small. The resulting moments if extrapolated to z 2 = 0, using a constant extrapolation at low z 2 , are in agreement with the results obtained by QCDSF within the errors of both computations.
Clearly, higher order perturbative matching is needed in order to be able to better estimate the systematic error of the perturbative correction. In addition computations at smaller lattice spacings are required in order to control the lattice spacing errors that should be enhanced for the shortest distance points. Finally, computations at larger volumes will provide a finer grid of ν values allowing for a more robust extraction of the derivatives of M(ν, z 2 ). This coupled with better statistics may allow us to extract higher moments as well. Detailed comparison between moments obtained from local and non-local matrix elements would require a control of all the above systematics as well as computations on the same ensembles. This is beyond the scope of the current work and will be done in the future.

IV. CONCLUSION
In this work we show that because the ratio M(ν, z 2 ) introduced in [16] is free of UV divergences and has a well defined continuum limit, it can be expanded into moments of parton distribution functions using OPE without any complications arising from power divergences due to the lattice regulator used to compute it non-perturbatively. In particular, if one expands this matrix element in lattice regularized twist-2 operators, which have power divergences due to breaking of the rotational symmetry on the lattice, the accompanying Wilson coefficients will also have power divergences which exactly cancel the power divergences of the matrix elements. Therefore, the re-summed OPE expansion of the reduced Ioffe time PDF is finite (as expected) due to cancellation of these power divergences order by order in the OPE expansion. Furthermore, using the relation between moments of PDFs and the derivatives of M(ν, z 2 ) with respect to the Ioffe time ν at ν = 0, it is shown that the moments of PDFs can be computed numerically from M(ν, z 2 ). The first two moments are found to agree within errors with those computed on the lattice by direct computations of matrix elements in the quenched approximation.
In the future, our calculations with dynamical quarks will improve on the systematics of the extraction of moments by addressing the sources of systematic errors that arise from the perturbative matching, the lattice spacing effects, as well as the numerical estimation of the derivatives of M(ν, z 2 ).