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 MS¯\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \overline{MS} $$\end{document} 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 agreement with the moments obtained from direct computations of local twist-2 matrix elements in the quenched approximation.


Introduction
It has been shown by Ji [1] that a numerical computation of parton distribution functions (PDFs) in Euclidean space using lattice QCD may be possible. Shortly after the basic idea was introduced several works [2][3][4][5][6] explored the properties of the methodology as well as introduced alternative approaches [7][8][9][10]. We refer the reader to [11] for a detailed review on the topic. However, in recent notes Rossi and Testa [12,13] raised a serious question about the validity of such approaches pointing out that the moments obtained from nonlocal operators are divergent. As a result these approaches may be impractical if one does not understand how to subtract non-perturbatively these unwanted effects. Despite the rebuttal in [14], in [13] it is argued that power divergences present in the moments of Ji's quasi-PDFs are a major obstruction in extracting the full PDF from lattice QCD calculations using this methodology.
In this work we will discuss potential power divergences in lattice QCD computations of Ioffe time parton distribution functions (PDFs) using the methodology introduced in [15,16]. Ioffe time PDFs are just the Fourier transforms of the longitudinal momentum fraction x PDFs, where the Ioffe time ν is the Fourier dual of the momentum fraction x. Ioffe time PDFs are directly related to the matrix elements computed in lattice QCD and therefore analyzing their behavior at finite lattice spacing is simpler. Furthermore, in the approach of [15,16], the light cone limit is taken as z 2 → 0 at fixed Ioffe time ν. On the other hand, the moments of PDFs are related to the coefficients of the Taylor expansion of an ultraviolet (UV) finite matrix element at zero Ioffe time and fixed z 2 . We use this property to compute the first two moments of PDFs in the M S scheme and find agreement with direct computations through local matrix elements of twist-2 operators. The issue of power divergent mixing of high dimensional operators with lower dimensional ones (the "infamous" trace operators) is a problem (actually an obstruction) if one wants to directly extract from lattice QCD calculations the PDFs as Fourier Transforms of hadronic matrix elements of the bi-local operator. On the contrary if, as it is done in the present paper, lattice data at short-distance are used to fit the operator product expansion (OPE) of the

JHEP11(2018)178
Ioffe time function in order to extract moments of PDFs, no question about power divergent mixing arises. A similar issue though in a different context was pointed out in [17]. In section 2, we remind the reader the details of the formalism, we discuss UV divergences and the expansion in moments. In section 3 we perform the computation of moments. section 4 contains our concluding remarks.

Ioffe time and pseudo parton distribution functions
Parton distribution functions can be computed from the hadronic matrix element of a quark bilinear operator with the quark and anti-quark fields separated by a finite distance. In the case of unpolarized, non-singlet parton densities the appropriate matrix element is where z is an arbitrary separation between the quark fields, p is an arbitrary momentum for the hadron,Ê(z, 0; A) is the z → 0 straight Wilson line in the fundamental representation, τ 3 is the flavor Pauli matrix, and γ α is a gamma matrix acting in spin space. Lorentz invariance dictates that this matrix element can be decomposed as 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) in light cone coordinates 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 to compute the PDF. 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 scale dependent PDFs [18][19][20]. 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 [21,22]. A convenient choice is z = (0, 0, 0, z 3 ) in Cartesian coordinates, α 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. 1 Considering a time local matrix element as above also allows for its computation in Euclidean space using lattice QCD. As

JHEP11(2018)178
it is well known [23], 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 [19,20,24,25]. 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 to consider 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,26], 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 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,26], the polynomial corrections for the case of the reduced function are smaller due to cancellation of certain polynomial terms in the numerator and the denominator. In particular, since the reduced function by construction satisfies M(0, z 2 ) = 1, one can see that B k (0) = 0 and In the following subsections we proceed with further discussion of UV divergences and consider the expansion in moments of the reduced function M(ν, z 2 ) from where the M S moments of PDFs can be extracted.

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 [27]). In this case, the renormalized matrix element JHEP11(2018)178 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 (2.8) 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 [28]) 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 (e.g. using the RI-MOM scheme [29]), and one 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.
In the limit of small space-like z 2 , the matrix elements of the numerator and the denominator of the ratio that defines M(ν, z 2 ) 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 in a scheme of our choice. The matrix elements can be computed in any regularization and renormalization scheme we desire provided that is matched to the perturbative scheme used to compute the Wilson coeficients. If we ignore higher twist effects, then the small z 2 expansion of our matrix element would read as

JHEP11(2018)178
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 ).
In conclusion the ratio M(ν, z 2 ) is a regularization and renormalization scheme independent quantity that can be expanded using OPE in the limit of small z 2 . Part of the polynomial corrections to the leading contribution is canceled by the polynomial terms in z 2 arising from the small z 2 expansion of the denominator. The above discussion can also be done using a cut-off scheme as a regulator. In this case power divergences will arise in both the computation of the Wilson coefficients and the computation of the matrix elements. However, divergences in the Wilson coefficients will cancel those of the matrix elements as M(ν, z 2 ) is regulator independent. This fact was pointed out in an other context in [30].

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 in the M S scheme in this section. First one can further simplify the expression in eq. (2.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. (2.12) we obtain where the product p 3 z 3 has been replaced by the Ioffe time ν. This formula for the moments is derived by the traditional definition a n (µ) =

JHEP11(2018)178
we can derive that where a n (µ), is the n-th moment of the parton distribution function. From this expression and eq. (2.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 from eq. (3.2) one can introduce m n as, Furthermore, eq. (3.5) implies that the Wilson coefficients are Since C(α, z 2 µ 2 , α s (µ)) is known analytically [20,25,31] 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. (3.7). 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 (3.10) In the above equations [· · · ] + denotes the "plus prescription". Integrating eq. (3.7) one obtains which are the well known leading order moments of the Altarelli-Parisi kernel, and d n = 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 (3.14) 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.
To illustrate this procedure we take as an example a recent quenched QCD calculation [16] which can be compared with the results from [32] where the moments were 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 computation is 0.093 fm. The same quenched theory was also used in [32] to study the moments of PDFs from direct computations of the corresponding twist-2 nucleon matrix elements, however in this case unimproved 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 figure 1, we plot the left hand side of eq. (3.6). 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

JHEP11(2018)178
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 figure 1, we plot the lowest two moments of the unpolarized isovector PDF computed using eq. (3.14) at scale µ = 3 GeV. These moments are plotted as function of z 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 comparable with the results obtained by QCDSF within the errors of both computations. It should be noted that the QCDSF results where computed at scale µ = √ 2 GeV, and for that reason we performed 2-loop running [33] of their results to the scale µ = 3 GeV where our moments were evaluated. In addition, the QCDSF computation was performed with unimproved Wilson fermions and therefore has O(a) errors. The remaining O(10%) difference between our results and those of QCDSF is expected within the systematic errors associated with perturbative matching, discretization errors, and the systematics of non-perturbative renomalization performed for the local twist-2 matrix elements.
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.

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 accompa-

JHEP11(2018)178
nying 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 be in agreement with those computed on the lattice by direct computations of matrix elements in the quenched approximation, within the statistical and systematic errors of the two calculations. Given that lattice calculations are much easier in the region of small Ioffe time ν, the methodology we presented here, which focuses on the small ν region, can lead to reliable non-perturbative computation of higher moments of PDFs. At this point we also wish to stress once again an important point. Namely, the issue of power divergent mixing of high dimensional operators with lower dimensional ones is a problem if one wants to directly extract from lattice QCD simulations the PDFs as Fourier transforms of hadronic matrix elements of the bilocal matrix element [34]. However, with the method employed in this study, where lattice data at short-distance are used to fit the OPE of the Ioffe time function in order to extract moments of PDFs, no issues with power divergent mixing arise.
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 ). In addition, computations on larger volumes will allow for better probing of the small Ioffe time ν. Such computations are currently being pursued and will be presented in future publications.

JHEP11(2018)178
Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract # DE-AC02-05CH11231.
Note added. While this work was in its completing stages, a preprint by A. Radyushkin [34] appeared on arXiv, where some of the points raised in this manuscript are also discussed.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.