# Moments of parton distribution functions of any order from lattice QCD

Andrea Shindler<EMAIL_ADDRESS>Institute for Theoretical
Particle Physics and Cosmology, TTK, RWTH Aachen University  Nuclear Science
Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
Department of Physics, University of California, Berkeley, CA 94720, USA

###### Abstract

We describe a procedure to determine moments of parton distribution functions
of any order in lattice QCD. The procedure is based on the gradient flow for
fermion and gauge fields. The flowed matrix elements of twist-2 operators
renormalize multiplicatively, and the matching with the physical matrix
elements can be obtained using continuum symmetries and the irreducible
representations of Euclidean $4$-dimensional rotations. We calculate the
matching coefficients at one-loop in perturbation theory for moments of any
order in the flavor non-singlet case. We also give specific examples of
operators that could be used in lattice QCD computations. It turns out that it
is possible to choose operators with identical Lorentz indices and still have
a multiplicative matching. One can thus use twist-2 operators exclusively with
temporal indices, thus substantially improving the signal-to-noise ratio in
the computation of the hadronic matrix elements.

††preprint: TTK-23-31

Introduction \- Quantum Chromodynamics (QCD), the theory of strong
interactions, comprehensively describes the interactions between quarks and
gluons across a broad energy spectrum, including the binding mechanisms
between these elementary particles. QCD is non-perturbative at low energies,
presenting a challenge in applying perturbative methods for the analysis of
hadron structures.

The study of hadron structure relies on the parton picture, where the dynamics
of quarks and gluons are characterized by parton distribution functions (PDF).
The precision of QCD calculations is intricately tied to the accuracy of these
PDFs, and this connection plays a crucial role in current measurements at the
LHC. It will continue to be essential for future investigations at the EIC,
HL-LHC, and LHeC. The current and future search for new physics depends
heavily on the precision of the comparison between experimental data and
theoretical predictions. This underscores the importance of a precise
determination of PDFs.

The understanding of collider processes based on QCD factorization allows for
the study of hadron structure using PDFs. Going beyond PDFs, generalized
parton distributions (GPDs), and transverse-momentum-dependent distributions
(TMDs) complement the 3-D description of hadrons. Although all the
distribution functions are not directly measurable, they can be determined by
factorizing the hard and soft parts of the cross-section of the process of
interest.

In deep inelastic scattering (DIS), a key role is played by the structure
functions, $F_{k}(x,Q^{2})~{}(k=1,2,3)$, which depend on $Q^{2}=-q^{2}$, where
$q$ is the 4-momentum transferred to the nucleon, and $x=Q^{2}/2(pq)$
represents the fraction of nucleon momentum carried by the parton. Here, $p$
denotes the nucleon momentum. For $Q^{2}$ much larger than the nucleon mass,
the factorization theorem [1]

$F_{k}(x,Q^{2})=\sum_{a}C_{k}^{a}(Q^{2}/\mu^{2})\otimes f_{a}(x,\mu^{2})\,,$
(1)

relates the structure functions to the convolution of the scale-dependent
parton distribution functions $f_{a}(x,\mu^{2})$ with the Wilson coefficients
$C_{k}^{a}(Q^{2}/\mu^{2})$ (summed over quark flavors and gluons, $a$). The
Wilson coefficients are calculable in perturbation theory, and the evolution
of the PDFs in $\mu$ is governed by the Dokshitzer-Gribov-Lipatov-Altarelli-
Parisi (DGLAP) equations [2, 3, 4, 5]. The DGLAP equations can be solved once
we specify an initial condition at a certain scale $\mu=\mu_{0}$. The
evolution of the PDFs is calculable in perturbation theory, but the precision
of experimental data always necessitates higher-order calculations. The
standard choice for renormalization and factorization scheme is the
$\overline{\text{MS}}$ scheme. PDFs can be determined experimentally by
analyzing a set of hard-scattering processes involving nucleons. There is a
long-standing community effort to determine PDFs using well-defined
parametrizations for their $x$-dependence and starting the evolution around
$\mu_{0}^{2}=1-4$GeV.

Lattice QCD can offer a theoretical input for the determination of the PDFs
and their evolution. The direct calculation of PDFs using lattice QCD poses
particular challenges due to the Euclidean geometry and the light-cone
dominance of the kinematics. The connection between PDFs and hadronic matrix
elements, which are calculable in lattice QCD, is established through the
moments of the PDFs, denoted as $\left\langle x^{n}\right\rangle$ [6, 7].
Lattice QCD calculations of the moments of the PDFs, pioneered in Refs. [8, 9,
10], provide, in principle, a means for the complete reconstruction of the
PDFs. This possibility has remained impractical due to the theoretical and
numerical challenges associated with computing high moments. Alternative
approaches have been developed to determine the $x$-dependence of the PDFs,
including the heavy-quarks OPE [11, 12], the quasi-PDF [13], the psuedo-PDF
[14], the OPE-based method of [15], the current-current approach [16], and the
hadron tensor method of [17].

Different high-energy scatterings are associated with different PDFs:
unpolarized, helicity, and transversity. Certain kinematical regions, as well
as specific PDFs, pose challenges for extraction from experimental data. This
emphasizes the significance of reconstructing the PDFs from lattice QCD
simulation data.

In this work, we revisit the possibility of calculating moments of PDFs,
which, if successful, would provide an important and complementary approach to
determine distribution functions from QCD. We describe a method that addresses
both the theoretical and numerical challenges faced in the past, which
hindered the calculation of moments of any order from lattice QCD.

Twist-2 operators \- QCD factorization, as expressed in Eq. (1), relates the
moments of the structure functions and the moments of the parton distribution
functions. Taking the example of $F_{1}$, the relation is given by e.g.

$\int_{0}^{1}dx~{}x^{n-1}F_{1}(x,Q^{2})=\sum_{a}C_{1,n}^{a}(Q^{2},\mu^{2})A_{n}^{a/h}(\mu)$
(2)

where $A_{n}^{a/h}(\mu)$ represents the reduced matrix element of certain
local operators, denoted by $O_{n}$, between the hadron states, denoted as
$h$, and $C_{1,n}^{a}(Q^{2},\mu^{2})$ are the corresponding Wilson
coefficients, calculable in perturbation theory. In the context of DIS
kinematics, the leading contribution, given in Eq. (2), arises from local
operators with the lowest twist, $\tau=d_{O}-n=2$, where $d_{O}$ is the
physical dimension of the local operator and $n$ is its spin. In the parton
model the reduced matrix element is related to the moments of the PDFs,
$A_{n}^{a/h}=\left\langle x^{n-1}\right\rangle_{a/h}$ [6, 7].

The form of the local operators $O_{n}$ can be constrained by requiring that
they belong to irreducible representations of the Lorentz group. This
guarantees that the renormalization procedure does not generate mixing with
operators of the same or lower dimensions. Lattice QCD simulations are
performed in a Euclidean setup though, and the symmetry group in a Euclidean
lattice is the hyper-cubic group H($4$). The reduced symmetry does not protect
the operators from a complicated mixing pattern, that gets more complicated
for larger $n$. Once the continuum limit has been performed one recovers
O($4$) symmetry.

We now briefly review how to construct irreducible representations of O($4$).
Although this is standard material (see, for example, Ref. [18]), it is useful
because it provides the basis to construct the operators that should be used
in lattice QCD computations.

A rank-$n$ tensor, denoted as $T_{\mu_{1}\cdots\mu_{n}}$, belongs to an
irreducible representation of the general linear group GL($4$) once the vector
indices are symmetrized

$T_{\left\\{\mu_{1}\cdots\mu_{n}\right\\}}=\frac{1}{n!}\sum_{\begin{subarray}{c}\sigma\in\text{all}\\\
\text{permutations}\end{subarray}}T_{\mu_{\sigma(1)}\cdots\mu_{\sigma(n)}}\,.$
(3)

In the case of the O($4$) group symmetry, tensors with given symmetries do not
represent irreducible representations because an additional operation is
allowed that commutes with the orthogonal transformations. If we contract $2$
indices, e.g. $(12)$, of the rank-$n$ tensor

$T^{(12)}_{\mu_{1}\cdots\mu_{n}}=T_{\alpha\alpha\mu_{3}\cdots\mu_{n}}=\delta_{\mu_{1}\mu_{2}}T_{\mu_{1}\cdots\mu_{n}}\,,$
(4)

we obtain a tensor with rank $n-2$ and the O($4$) transformation commutes with
the operation of contraction. We can repeat the contraction operation for any
pair of vector indices, obtaining $n(n-1)/2$ traces. The subspace of the
rank-$n$ tensors for which all pairs of traces are zero is invariant under
O($4$) rotation, and it can be shown that every rank-$n$ tensor can always be
decomposed in a tensor with vanishing trace, denoted as $\widehat{T}$, and a
linear combination of $n(n-1)/2$ terms of all possible traces

$T_{\mu_{1}\cdots\mu_{n}}=\widehat{T}_{\mu_{1}\cdots\mu_{n}}+\delta_{\mu_{1}\mu_{2}}T^{(12)}_{\mu_{1}\cdots\mu_{n}}+\cdots\,.$
(5)

This decomposition is invariant under O($4$) transformations. For $n=2$ the
decomposition is very easy to obtain and the traceless tensor reads

$\widehat{T}_{\mu_{1}\mu_{2}}=T_{\mu_{1}\mu_{2}}-\frac{1}{4}\delta_{\mu_{1}\mu_{2}}T_{\alpha\alpha}\,.$
(6)

For $n=3$, we have $3$ different traces, and the traceless part of the tensor
can be computed in the same way

$\displaystyle\widehat{T}_{\mu_{1}\mu_{2}\mu_{3}}=T_{\mu_{1}\mu_{2}\mu_{3}}-$
(7) $\displaystyle-$
$\displaystyle\frac{1}{6}\left[\delta_{\mu_{1}\mu_{2}}T_{\alpha\alpha\mu_{3}}+\delta_{\mu_{1}\mu_{3}}T_{\alpha\mu_{2}\alpha}+\delta_{\mu_{2}\mu_{3}}T_{\mu_{1}\alpha\alpha}\right]\,.$

Traceless tensors are invariant under vector index permutations and can be
used as a starting point to construct all the irreducible representations of
O($4$) using Young symmetrizers [18]. For the purpose of this work, the only
relevant result is that traceless and symmetrized rank-$n$ tensors are an
irreducible representation of O($4$).

We restrict ourselves to unpolarized structure functions, meaning that we
consider hadronic matrix elements averaged over the hadron polarization. Other
moments of different distributions can be studied using the same procedure
described here. We further limit our consideration to flavor non-singlet
operators to avoid mixing with gluonic operators, but the same method is
applicable to singlet operators.

We consider the following set of operators

$O_{n}^{rs}(x)=O^{rs}_{\mu_{1}\cdots\mu_{n}}(x)=\overline{\psi}^{r}(x)\gamma_{\left\\{\mu_{1}\right.}\overset{\leftrightarrow}{D}_{\mu_{2}}\cdots\overset{\leftrightarrow}{D}_{\left.\mu_{n}\right\\}}\psi^{s}(x)\,.$
(8)

with quark fields, $\psi^{r}$, of different flavors, $r\neq s$, and the
covariant derivative $D_{\mu}=\partial_{\mu}+G_{\mu}$, where $G_{\mu}$
represent the gluon field. The symmetrization of the covariant derivative
$\overset{\leftrightarrow}{D}_{\mu}=\frac{1}{2}\left({\overset{\rightarrow}{D}_{\mu}}-{\overset{\leftarrow}{D}_{\mu}}\right)$
guarantees a definite transformation under charge conjugation.

In the continuum, traceless non-singlet twist-2 operators, with symmetrized
Lorentz indices, renormalize multiplicatively. For example, in the MS scheme,
the bare operator $O_{n,{\text{B}}}^{rs}$ renormalizes as follows

$O_{n}^{rs}(x)=Z^{\text{MS}}_{n}O_{n,{\text{B}}}^{rs}(x)\,,$ (9)

with renormalization constant $Z^{\text{MS}}_{n}$. The intrinsic non-
perturbative nature of the hadronic matrix elements suggests the use of
lattice QCD. However, the breaking of the rotational group symmetry into the
hypercubic group H($4$) adds complications to the renormalization.

On an hypercubic lattice, the Lorentz indices of the twist-2 operators must be
chosen while considering the reduced hypercubic symmetry group H($4$).
Irreducible representations of O($4$) generally become reducible
representations of H($4$) inducing unwanted mixings under renormalization [8,
9, 19, 20]. The irreducible representations of H($4$) allow mixings with
lower-dimensional operators, and, additionally, complicate the mixing with
operators of the same dimensions. A notable example of this complication
occurs with the operator $O_{3}$. In the continuum, after symmetrizing over
the Lorentz indices and subtracting the traces to obtain $\widehat{O}_{3}$,
the renormalization is multiplicative (e.g. in the MS scheme). However, on the
lattice, terms like $1/a^{2}\delta_{\mu_{i}\mu_{j}}\cos(ap_{\mu_{j}})$ are
allowed by H($4$) symmetry and cannot be removed even after subtracting the
trace [8]. This pertains to the fact that $O_{3}$ belongs to different
irreducible representations of H($4$) depending on the choice of the Lorentz
indices. The operator $O_{3}$ with all the indices equal, i.e.
$\mu_{1}=\mu_{2}=\mu_{3}$, mixes with the lower-dimensional vector current,
causing a power divergence proportional to $1/a^{2}$. For the operators with
$\mu_{1}\neq\mu_{2}=\mu_{3}$, such as the irreducible representation
$O_{411}-O_{433}$, the renormalization is multiplicative.111 Strictly
speaking, $O_{411}-O_{433}$ belongs to an irreducible representation of H($4$)
that appears more than once. Therefore, a mixing between operators belonging
to these equivalent representations can occur [19]. In perturbation theory the
mixing seems to be numerically small.  We note that while $O_{411}$ is also
affected by power divergences, the subtraction $O_{411}-O_{433}$ guarantees
their removal. Another choice is $\mu_{1}\neq\mu_{2}\neq\mu_{3}$. In this
case, the operator belongs to another irreducible representation of H($4$) and
renormalizes multiplicatively without introducing power divergences. However,
this last example, while optimal in terms of renormalization, is far from
optimal when calculating hadronic matrix elements222In Eq. (10) for simplicity
we omit the flavor indices.

$\left\langle
h(p)|O_{\left\\{\mu_{1}\cdots\mu_{n}\right\\}}|h(p)\right\rangle\
=p_{\mu_{1}}\cdots
p_{\mu_{n}}A_{n}^{h}(\mu)~{}+~{}\delta_{\mu_{i}\mu_{j}}{\text{terms}}\,,$ (10)

because it requires non-vanishing external momentum in at least $2$ spatial
directions. This substantially degrades the quality of the signal-to-noise
ratio that is realized in numerical lattice QCD simulations. For
$O_{411}-O_{433}$, one needs a non-vanishing external momentum in at least $1$
spatial directions, and despite the subtraction of power divergences, one
might hope to achieve a better signal-to-noise ratio. The operator with all
the same indices has a power divergence that requires a complicated non-
perturbative procedure to be subtracted and is not a viable solution. This is
unfortunate because the matrix element of $O_{444}$, for example, does not
require any external spatial momentum and should be optimal in terms of the
signal-to-noise ratio.

The mixing pattern becomes more cumbersome for higher dimensional operators,
$n\geq 4$. For $n=4$, choosing all $4$ indices differently results in the only
$2$ irreducible representations: the totally symmetric and anti-symmetric
representations. One then needs a non-vanishing external momentum in all $3$
spatial directions to calculate the hadronic matrix elements[19]. For $n>4$,
mixings are unavoidable. The complicated renormalization pattern just
described, together with the need of an external momentum to calculate the
hadronic matrix elements, has, de facto, prevented the lattice calculation of
the moments of parton distribution functions for $n>4$ and rendered the
calculations for $n=3,4$ challenging.

A very interesting idea to circumvent the renormalization issue, based on the
recovery of the continuum rotational symmetry, was proposed in Ref. [21]. Even
though the method was not further pursued numerically, it represents a first
important step toward resolving the theoretical challenges encountered when
calculating moments of PDFs.

Gradient Flow \- The gradient flow (GF) for gauge and fermion fields [22, 23,
24] provides a different way to regulate short-distance singularities of the
twist-2 operators. The connection with the physical renormalized matrix
elements at vanishing flow time, $t=0$, is obtained with a short-flow-time
expansion (SFTX) after performing the continuum limit of the hadronic matrix
elements at a fixed physical value of the flow time $t$ [25].

The flowed twist-2 operators, written in terms of flowed fermions,
$\chi^{r}(x,t)$ and $\overline{\chi}^{r}(x,t)$, and gauge fields,
$B_{\mu}(x,t)$, are given by

$O_{n}^{rs}(x,t)=\overline{\chi}^{r}(x,t)\gamma_{\left\\{\mu_{1}\right.}\overset{\leftrightarrow}{D}_{\mu_{2}}\cdots\overset{\leftrightarrow}{D}_{\left.\mu_{n}\right\\}}\chi^{s}(x,t)\,,$
(11)

where $\left\\{\cdots\right\\}$ denotes symmetrization (cf. Eq. (3)). These
operators are particularly advantageous because they renormalize
multiplicatively also on the lattice with a renormalization factor that
depends only on the fermion content of the operator,

$O_{n}^{rs}(t)=Z_{n}O_{n,B}^{rs}(t)\,,\quad Z_{n}=Z_{\chi}\,,$ (12)

where $Z_{\chi}^{1/2}$ is the renormalization constant of the flowed fermion
fields [24]. For correlation functions containing flowed twist-2 operators,
the flow-time $t$ provides a regulator for short-distance singularities.
Beside the renormalization of the bare parameters of the theory, the flowed
fermion fields, and any local field at $t=0$, correlation functions do not
require any additional renormalization. For lattice QCD applications, it is
convenient to define $Z_{\chi}$ in a scheme that is regularization
independent, as then the matching and the renormalization can be performed in
the same scheme. The choice of the specific scheme is not relevant for the
general discussion, but it becomes important when calculating the matching
coefficients. In the calculations of the matching coefficients described
below, we renormalize the twist-2 operators using the so-called ringed fields
[26], defined by the SU($N_{c}$) gauge invariant and regularization
independent condition

$\left\langle\mathring{\overline{\chi}}_{r}(x,t){\overset{\leftrightarrow}{\not{D}}}\mathring{\chi}_{r}(x,t)\right\rangle=-\frac{N_{c}}{(4\pi)^{2}t^{2}}\,.$
(13)

Imposing this condition in dimensional regularization, with $D=4-2\epsilon$,
leads to a finite renormalization between the ringed fields and renormalized
fields in the MS scheme

$\displaystyle\chi_{r}(x,t)$ $\displaystyle=$ $\displaystyle\left(8\pi
t\right)^{\epsilon/2}\zeta_{\chi}^{1/2}\mathring{\chi}_{r}(x,t)$
$\displaystyle\overline{\chi}_{r}(x,t)$ $\displaystyle=$
$\displaystyle\left(8\pi
t\right)^{\epsilon/2}\zeta_{\chi}^{1/2}\mathring{\overline{\chi}}_{r}(x,t)\,.$
(14)

The 1-loop result for the finite renormalization is given by [26]

$\zeta_{\chi}=1-\frac{\overline{g}^{2}}{(4\pi)^{2}}C_{F}\left(3\log(8\pi\mu^{2}t)-\log(432)\right)\,,$
(15)

where $\overline{g}$ is the strong coupling renormalized in the MS scheme at
the renormalization scale $\mu$, and $C_{F}=N_{c}^{2}-1/2N_{c}$. This result
has been extended to O($\overline{g}^{4}$) in Refs. [27, 28]. Results obtained
in the MS scheme can be converted to the $\overline{\text{MS}}$ scheme by
replacing the renormalization scale $\mu$ with

$\mu^{2}=\overline{\mu}^{2}\frac{{\rm{e}}^{\gamma_{E}}}{(4\pi)}\,,$ (16)

where $\overline{\mu}$ is the renormalization scale in the
$\overline{\text{MS}}$ scheme.

O($a$) improvement \- The lattice QCD calculation of hadronic matrix elements
of twist-2 operators is affected by cutoff effects. For non-perturbatively
improved clover fermions, the hadronic matrix elements of the standard twist-2
operators (with $t=0$) are still affected by O($a$) cutoff effects that could
be removed by adopting an O($a$) improved definition of the twist-2 operators.
Excluding the case $n=2$, where the corresponding improvement coefficients
have been determined at 1-loop in perturbation theory [29], the continuum
limit is reached with O($a$) discretization effects.

If we consider flowed twist-2 operators for a generic $n$, the hadronic matrix
elements are affected only by O($am$) cutoff effects, where $m$ is the quark
mass. Additional improvement coefficients specific to the twist-2 operators
are not needed. In the discussion of the method’s applications below, we also
consider ratios of flowed matrix elements where both $Z_{\chi}$ and the
O($am$) cutoff effects are eliminated. In summary, regardless of the lattice
setup adopted, the use of flowed twist-2 operators provides a significant
simplification in the renormalization pattern. Moreover, for non-
perturbatively improved clover fermions, it also improves the parametric
scaling towards the continuum limit.

Matching coefficients \- After performing the continuum limit of the hadronic
matrix element for a fixed $t>0$ value of the flow time, the connection with
physical renormalized matrix element at $t=0$ is obtained with a SFTX [25].
The SFTX of the twist-2 operators $O_{n}(x,t)$, in general contains power
divergences with terms like $1/t^{m}\delta_{\mu_{i}\mu_{j}}$, where
$m=(d_{n}-d_{p})/2$ depends on the dimensions of $O_{n}$, $d_{n}$, and on the
dimension of the lower-dimensional operator, $d_{p}$. These terms, though, are
classified using the continuum O($4$) symmetry and, as a result, the SFTX of
the symmetrized and traceless operators, $\widehat{O}_{n}$, receives
contributions only from the corresponding traceless operators at $t=0$.

We renormalize the flowed fermion fields adopting ringed fields, introduced in
the previous section, and define the ringed twist-2 operators

$\displaystyle\widehat{\mathring{O}}_{n}^{rs}(x,t)$ $\displaystyle=$
$\displaystyle\mathring{\overline{\chi}}^{r}(x,t)\gamma_{\left\\{\mu_{1}\right.}\overset{\leftrightarrow}{D}_{\mu_{2}}\cdots\overset{\leftrightarrow}{D}_{\left.\mu_{n}\right\\}}\mathring{\chi}^{s}(x,t)-$
(17) $\displaystyle-$
$\displaystyle\text{terms~{}with~{}}\delta_{\mu_{i}\mu_{j}}\,,$

where the subtraction guarantees that the resulting operator has vanishing
traces. The SFTX of $\widehat{\mathring{O}}_{n}^{rs}$ is

$\widehat{\mathring{O}}_{n}^{rs}(t)=c_{n}(t,\mu)\widehat{O}_{n,\text{MS}}^{rs}(\mu)+\cdots\,,$
(18)

where we make explicit only the dependence on the flow time and the
renormalization scale. The neglected terms are higher-order contributions of
positive powers of $t$.

The matching coefficients are determined considering off-shell amputated one-
particle irreducible (1PI) Green’s functions containing the flowed operators
$\widehat{\mathring{O}}_{n}^{rs}(t)$. The matching equations with massless
fermions

$\left\langle\psi^{r}\widehat{\mathring{O}}_{n}^{rs}(t)\overline{\psi}^{s}\right\rangle=c_{n}(t,\mu)\left\langle\psi^{r}\widehat{O}_{n,\text{MS}}^{rs}(t=0,\mu)\overline{\psi}^{s}\right\rangle$
(19)

are solved in $D=4-2\epsilon$ at 1-loop in perturbation theory with

$c_{n}(t,\mu)=1+\frac{\overline{g}^{2}(\mu)}{\left(4\pi\right)^{2}}c_{n}^{(1)}(t,\mu)+O(\overline{g}^{4})\,.$
(20)

The matching coefficients $c_{n}(t,\mu)$ depend on the renormalized coupling
and are independent of the soft scales. For their calculation, one uses
standard techniques (see for example Ref. [30] and Refs. therein), expanding
the integrands of the loop integrals of the matching equations (19) in all
soft scales keeping $t$ fixed. The $t=0$ 1-loop contribution on the r.h.s of
Eq. (19) vanishes because the expansion leads to scaleless, i.e. vanishing,
integrals in D-dimensions. After renormalizing the gauge coupling and the
flowed fermion fields, the left-hand side is UV finite. The infrared
singularities stemming from the expansion in the soft scales, and regulated in
dimensional regularization, match exactly the UV poles of $Z_{n}^{\text{MS}}$.
The finite matching coefficients are the result of this procedure.

The Feynman diagrams contributing are depicted 333To draw the Feynman diagrams
we have used FeynGame [31]. in Fig. 1. The cross indicates the insertion of
the twist-2 operators, the single straight and wavy lines respectively
represent the fermion and gluon propagator, and the double straight lines
represent the kernel of the fermionic GF operator. The open circle represents
a GF vertex, and the rest are standard QCD vertices. The Feynman rules for the
operators are the same as in the unflowed case, and the remaining Feynman
rules related to the expansion of the GF equation can be found, for example,
in Refs. [32, 30].

Figure 1: Feynman diagram for the O($\overline{g}^{2}$) calculations of the
matching coefficients in Eq. (19). The circle with the cross represents the
twist-2 operator and the empty circles the GF vertices. The double line
represents the fermionic kernel of the GF differential operator. The arrows
beside the fermion lines and the kernel lines indicate the direction of
increasing flow time.

The final result for the matching coefficients is

$c_{n}^{(1)}(t,\mu)=C_{F}\left[\gamma_{n}\log\left(8\pi\mu^{2}t\right)+B_{n}\right]\,,$
(21)

where

$\gamma_{n}=1+4\sum_{j=2}^{n}\frac{1}{j}-\frac{2}{n(n+1)}\,,$ (22)

and

$\displaystyle B_{n}$ $\displaystyle=$
$\displaystyle\frac{4}{n(n+1)}+4\frac{n-1}{n}\log
2+\frac{2-4n^{2}}{n(n+1)}\gamma_{E}-$ (23) $\displaystyle-$
$\displaystyle\frac{2}{n(n+1)}\psi(n+2)+\frac{4}{n}\psi(n+1)-4\psi(2)-$
$\displaystyle-$ $\displaystyle
4\sum_{j=2}^{n}\frac{1}{j(j-1)}\frac{1}{2^{j}}\phi(1/2,1,j)-\log\left(432\right)\,.$

In this expression $\gamma_{E}=0.57221\ldots$ is the Euler’s constant,
$\psi(z)$ is the digamma function, and $\phi(z,s,a)$ is the Lerch transcendent
defined as $\phi(z,s,a)=\sum_{k=0}^{\infty}\frac{z^{k}}{(k+a)^{s}}$.

The expression for $\gamma_{n}$ matches the 1-loop anomalous dimension of the
twist-2 operators [33] and provides a welcome check of the calculation. For
$n>2$ the finite part, $B_{n}$, in Eq. (23) provides a new result. The result
with $n=2$ was already obtained in Ref. [26] in the analysis of the energy-
momentum tensor, and we reproduce their result.

The fact that the matching of the flowed matrix elements is so simple opens
the possibility of computing moments of the parton distribution functions of
any order. A possible operational procedure is the following: After fixing
$n$, one constructs the symmetric and traceless operator
$\widehat{O}_{n}^{rs}$. One proceeds calculating the matrix elements between
hadron states of

$\left\langle h(p)|\widehat{\mathring{O}}_{n}^{rs}(t)|h(p)\right\rangle\,,$
(24)

using standard techniques based on the spectral decomposition. To apply the
decomposition in the lattice calculation of the 3-point function, one needs to
ensure that the physical distance between the interpolators of the hadron
states and of the flowed operators is much larger than the flow-time radius
$\sqrt{8t}$. After performing the continuum limit, the renormalized matrix
element in the MS scheme is calculated simply by

$A_{n}^{\text{MS}}(\mu)=c_{n}(t,\mu)^{-1}A_{n}(t)\,.$ (25)

This matching is correct only if the flowed fields are defined as ringed
fields. If the flowed fields are renormalized in a different scheme, the
finite part of the matching coefficient has to be modified.

If the matching is successful, the l.h.s of Eq. (25) should be independent of
the flow time. Violations stem from higher dimensional operators as linear
terms in $t$, and from higher order terms in the perturbative expansion of the
matching coefficients. To mitigate these possible systematics it is important
to extend to O($\overline{g}^{4}$) the calculation presented here and conduct
a thorough numerical study of the residual flow-time dependence of the
hadronic matrix element after the matching.

Applications \- The simplicity of the matching is guaranteed once the
operators are symmetrized on the Lorentz indices, and all the traces have been
subtracted. In the case of $n=3$, the traceless operator is

$\displaystyle\widehat{O}_{3}^{rs}=\widehat{O}_{\left\\{\mu_{1}\mu_{2}\mu_{3}\right\\}}^{rs}=O_{\left\\{\mu_{1}\mu_{2}\mu_{3}\right\\}}^{rs}-$
$\displaystyle-\frac{1}{6}\left[\delta_{\mu_{1}\mu_{2}}O^{rs}_{\left\\{\alpha\alpha\mu_{3}\right\\}}+\delta_{\mu_{1}\mu_{3}}O^{rs}_{\left\\{\alpha\mu_{2}\alpha\right\\}}+\delta_{\mu_{2}\mu_{3}}O^{rs}_{\left\\{\mu_{1}\alpha\alpha\right\\}}\right]\,,$
(26)

where repeated indices are summed over. This implies that for example we can
now choose $\mu_{1}=\mu_{2}=\mu_{3}=4$ and calculate the matrix element of

$\widehat{O}_{444}=\frac{1}{2}\left[O_{444}-\frac{1}{3}\left(O_{kk4}+O_{k4k}+O_{4kk}\right)\right]\,.$
(27)

The hadronic matrix element of $\widehat{O}_{444}(t)$ is calculable with
vanishing external momentum and, at finite flow time, it renormalizes
multiplicatively. The matching with the MS scheme at renormalization scale
$\mu$ is obtained using Eq. (25). In a similar way, traceless operators for
any $n$ can be constructed. To exemplify how to construct traceless operators,
for $n=4$ we have

$\widehat{O}_{4444}=O_{4444}-\frac{3}{4}O_{\left\\{\alpha\alpha
44\right\\}}+\frac{1}{16}O_{\left\\{\alpha\alpha\beta\beta\right\\}}\,.$ (28)

Again the hadronic matrix element of $\widehat{O}_{4444}(t)$ is calculable
with no external spatial momentum, and the matching with the MS operator is
multiplicative.

The continuum limit of the reduced hadronic matrix element of
$\widehat{\mathring{O}}_{n}$, denoted as $A_{n}^{h}(t)=\left\langle
x^{n-1}\right\rangle^{h}(t)$, would require a determination of $Z_{\chi}$ from
lattice QCD simulations. However, this can be avoided taking the second moment
$A_{2}^{h,\text{MS}}(\mu)=\left\langle x\right\rangle_{\text{MS}}^{h}(\mu)$ as
already computed with any method of choice. Then all the other moments for
$n>2$ are calculable from

$\left\langle x^{n-1}\right\rangle_{\text{MS}}^{h}(\mu)=\left\langle
x\right\rangle_{\text{MS}}^{h}(\mu)\frac{c_{2}(t,\mu)}{c_{n}(t,\mu)}\frac{\left\langle
x^{n-1}\right\rangle^{h}(t)}{\left\langle x\right\rangle^{h}(t)}\,,$ (29)

where $c_{2}(t,\mu)$ and $c_{n}(t,\mu)$ are given by Eqs. (21-23).

The ratios $\frac{\left\langle x^{n-1}\right\rangle^{h}(t)}{\left\langle
x\right\rangle^{h}(t)}$ have a finite continuum limit, and there is no need to
renormalize the flowed fields. Moreover, these ratios exhibit O($a^{2}$)
scaling violations once the lattice theory has been non-pertubatively
improved. If convenient numerically, it is also possible to use different
hadrons for the $n-$th and the second moment matrix elements in Eq. (29)
because the matching coefficients are independent on the hadron for which we
calculate the matrix element.

The gradient flow has been previously employed to define quasi-PDF in [34],
and the 1-loop matching has been calculated in [35]. While, in principle, the
determination of moments of PDFs should allow a precise reconstruction of the
PDFs, it would be intriguing to explore whether the determination of the PDF
moments outlined in this work can be integrated with quasi-PDF calculations.

Summary \- We have described a method that resolves both the theoretical and
numerical challenges faced when calculating moments of parton distribution
functions from lattice QCD. A lattice regulator breaks the Euclidean
rotational symmetry into the hypercubic group. As a result, the
renormalization of twist-2 operators on the lattice is a substantial hurdle to
determine moments of parton distribution functions $\left\langle
x^{n-1}\right\rangle$. For $n>4$, power divergences in the lattice spacing are
unavoidable, and for $n=3,4$ the cancellation of power divergences is still
possible by choosing suitable irreducible representations of the hypercubic
group. The price to pay is that the calculation of the hadronic matrix
elements necessitates the use of a spatial external momentum, worsening the
signal-to-noise ratio of the correlation functions. At finite flow time $t$,
the correlation functions of twist-2 operators are finite once the bare
parameters of the theory and the flowed fermion fields are renormalized. After
performing the continuum limit, the matching with the physical matrix element
can be done at short flow time using the continuum symmetries of the theory.
Using traceless and symmetrized operators guarantees that the matching is
multiplicative, and we have calculated, at 1-loop in perturbation theory, the
matching coefficients. This result opens the way to determine moments of any
order of the parton distribution functions. While, in principle, this enables
the reconstruction of PDFs from lattice QCD simulations, it would be very
interesting to investigate how the approach described in this work can
complement recent efforts aimed at the direct determination of PDFs from
lattice QCD correlators.

Acknowledgments I would like to thank Andre Walker-Loud for constant
encouragement and a critical reading of the manuscript. I also thank R.
Harlander and T. Luu for suggestions on how to improve the original draft. I
have benefited from many discussions about the gradient flow with R.
Harlander, J. Kim, T. Luu, E. Mereghetti, C. Monahan, M. Rizik, P. Stoffer and
A. Walker-Loud. I want to thank the Nuclear Theory Group of the Lawrence
Berkeley National Lab and the University of California Berkeley for their
hospitality. I acknowledge funding support from Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) through grant
513989149 and under the National Science Foundation grant PHY-2209185.

## References

  * Collins  _et al._ [1989] J. C. Collins, D. E. Soper, and G. F. Sterman, Factorization of Hard Processes in QCD, Adv. Ser. Direct. High Energy Phys. 5, 1 (1989), arXiv:hep-ph/0409313 .
  * Gribov and Lipatov [1972] V. N. Gribov and L. N. Lipatov, Deep inelastic e p scattering in perturbation theory, Sov. J. Nucl. Phys. 15, 438 (1972).
  * Lipatov [1974] L. N. Lipatov, The parton model and perturbation theory, Yad. Fiz. 20, 181 (1974).
  * Dokshitzer [1977] Y. L. Dokshitzer, Calculation of the Structure Functions for Deep Inelastic Scattering and e+ e- Annihilation by Perturbation Theory in Quantum Chromodynamics., Sov. Phys. JETP 46, 641 (1977).
  * Altarelli and Parisi [1977] G. Altarelli and G. Parisi, Asymptotic Freedom in Parton Language, Nucl. Phys. B 126, 298 (1977).
  * Curci  _et al._ [1980] G. Curci, W. Furmanski, and R. Petronzio, Evolution of Parton Densities Beyond Leading Order: The Nonsinglet Case, Nucl. Phys. B 175, 27 (1980).
  * Collins and Soper [1982] J. C. Collins and D. E. Soper, Parton Distribution and Decay Functions, Nucl. Phys. B 194, 445 (1982).
  * Kronfeld and Photiadis [1985] A. S. Kronfeld and D. M. Photiadis, Phenomenology on the Lattice: Composite Operators in Lattice Gauge Theory, Phys. Rev. D 31, 2939 (1985).
  * Martinelli and Sachrajda [1987] G. Martinelli and C. T. Sachrajda, Pion Structure Functions From Lattice QCD, Phys. Lett. B 196, 184 (1987).
  * Martinelli and Sachrajda [1988] G. Martinelli and C. T. Sachrajda, A Lattice Calculation of the Pion’s Form-Factor and Structure Function, Nucl. Phys. B 306, 865 (1988).
  * Aglietti  _et al._ [1998] U. Aglietti, M. Ciuchini, G. Corbo, E. Franco, G. Martinelli, and L. Silvestrini, Model independent determination of the shape function for inclusive B decays and of the structure functions in DIS, Phys. Lett. B 432, 411 (1998), arXiv:hep-ph/9804416 .
  * Detmold and Lin [2006] W. Detmold and C. J. D. Lin, Deep-inelastic scattering and the operator product expansion in lattice QCD, Phys. Rev. D 73, 014501 (2006), arXiv:hep-lat/0507007 .
  * Ji [2013] X. Ji, Parton Physics on a Euclidean Lattice, Phys. Rev. Lett. 110, 262002 (2013), arXiv:1305.1539 [hep-ph] .
  * Radyushkin [2017] A. Radyushkin, Nonperturbative Evolution of Parton Quasi-Distributions, Phys. Lett. B 767, 314 (2017), arXiv:1612.05170 [hep-ph] .
  * Chambers  _et al._ [2017] A. J. Chambers, R. Horsley, Y. Nakamura, H. Perlt, P. E. L. Rakow, G. Schierholz, A. Schiller, K. Somfleth, R. D. Young, and J. M. Zanotti, Nucleon Structure Functions from Operator Product Expansion on the Lattice, Phys. Rev. Lett. 118, 242001 (2017), arXiv:1703.01153 [hep-lat] .
  * Ma and Qiu [2018] Y.-Q. Ma and J.-W. Qiu, Extracting Parton Distribution Functions from Lattice QCD Calculations, Phys. Rev. D 98, 074021 (2018), arXiv:1404.6860 [hep-ph] .
  * Liang  _et al._ [2020] J. Liang, T. Draper, K.-F. Liu, A. Rothkopf, and Y.-B. Yang (XQCD), Towards the nucleon hadronic tensor from lattice QCD, Phys. Rev. D 101, 114503 (2020), arXiv:1906.05312 [hep-ph] .
  * Hamermesh [1962] M. Hamermesh,  _Group Theory and Its Application to Physical Problems_, Addison-Wesley series in physics (Addison-Wesley Publishing Company, 1962).
  * Beccarini  _et al._ [1995] G. Beccarini, M. Bianchi, S. Capitani, and G. Rossi, Deep inelastic scattering in improved lattice QCD. 2. The second moment of structure functions, Nucl. Phys. B 456, 271 (1995), arXiv:hep-lat/9506021 .
  * Gockeler  _et al._ [1996] M. Gockeler, R. Horsley, E.-M. Ilgenfritz, H. Perlt, P. E. L. Rakow, G. Schierholz, and A. Schiller, Lattice operators for moments of the structure functions and their transformation under the hypercubic group, Phys. Rev. D 54, 5705 (1996), arXiv:hep-lat/9602029 .
  * Davoudi and Savage [2012] Z. Davoudi and M. J. Savage, Restoration of Rotational Symmetry in the Continuum Limit of Lattice Field Theories, Phys. Rev. D 86, 054505 (2012), arXiv:1204.4146 [hep-lat] .
  * Lüscher, M. [2010] Lüscher, M., Properties and uses of the Wilson flow in lattice QCD, JHEP 1008, 071, arXiv:1006.4518 [hep-lat] .
  * Lüscher, M. and Weisz [2011] Lüscher, M. and P. Weisz, Perturbative analysis of the gradient flow in non-abelian gauge theories, JHEP 1102, 051, arXiv:1101.0963 [hep-th] .
  * Lüscher, M. [2013] Lüscher, M., Chiral symmetry and the Yang–Mills gradient flow, JHEP 1304, 123, arXiv:1302.5246 [hep-lat] .
  * Lüscher [2014] M. Lüscher, Future applications of the Yang-Mills gradient flow in lattice QCD, PoS LATTICE2013, 016 (2014), arXiv:1308.5598 [hep-lat] .
  * Makino and Suzuki [2014] H. Makino and H. Suzuki, Lattice energy–momentum tensor from the Yang-Mills gradient flow–inclusion of fermion fields, PTEP 2014, 063B02 (2014), arXiv:1403.4772 [hep-lat] .
  * Harlander  _et al._ [2018] R. V. Harlander, Y. Kluth, and F. Lange, The two-loop energy–momentum tensor within the gradient-flow formalism, Eur. Phys. J. C78, 944 (2018), [Erratum: Eur. Phys. J.C79,no.10,858(2019)], arXiv:1808.09837 [hep-lat] .
  * Artz  _et al._ [2019] J. Artz, R. V. Harlander, F. Lange, T. Neumann, and M. Prausa, Results and techniques for higher order calculations within the gradient-flow formalism, JHEP 06, 121, [erratum: JHEP10,032(2019)], arXiv:1905.00882 [hep-lat] .
  * Capitani  _et al._ [2001] S. Capitani, M. Gockeler, R. Horsley, H. Perlt, P. E. L. Rakow, G. Schierholz, and A. Schiller, Renormalization and off-shell improvement in lattice perturbation theory, Nucl. Phys. B 593, 183 (2001), arXiv:hep-lat/0007004 .
  * Mereghetti  _et al._ [2022] E. Mereghetti, C. J. Monahan, M. D. Rizik, A. Shindler, and P. Stoffer, One-loop matching for quark dipole operators in a gradient-flow scheme, JHEP 04, 050, arXiv:2111.11449 [hep-lat] .
  * Harlander  _et al._ [2020] R. V. Harlander, S. Y. Klein, and M. Lipp, FeynGame, Comput. Phys. Commun. 256, 107465 (2020), arXiv:2003.00896 [physics.ed-ph] .
  * Rizik  _et al._ [2020] M. D. Rizik, C. J. Monahan, and A. Shindler (SymLat), Short flow-time coefficients of $CP$-violating operators, Phys. Rev. D 102, 034509 (2020), arXiv:2005.04199 [hep-lat] .
  * Gross and Wilczek [1974] D. J. Gross and F. Wilczek, ASYMPTOTICALLY FREE GAUGE THEORIES. 2., Phys. Rev. D 9, 980 (1974).
  * Monahan and Orginos [2017] C. Monahan and K. Orginos, Quasi parton distributions and the gradient flow, JHEP 03, 116, arXiv:1612.01584 [hep-lat] .
  * Monahan [2018] C. Monahan, Smeared quasidistributions in perturbation theory, Phys. Rev. D97, 054507 (2018), arXiv:1710.04607 [hep-lat] .