  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.
