We present a new method to calculate moments of parton distribution functions of any order with
lattice QCD computations. This method leverages the gradient flow for fermion and gauge fields.
The flowed matrix elements of twist-2 operators renormalize multiplicatively, and the matching
with physical matrix elements is achieved through the use of continuum symmetries. We derive
the matching coefficients at one-loop in perturbation theory for moments of any order in the flavor
non-singlet case and provide specific examples of operators suitable for lattice QCD computations.
The multiplicative renormalization and matching are independent of the choice of Lorentz indices,
allowing the use of temporal indices for twist-2 operators of any dimension. This approach should
then also significantly enhance the signal-to-noise ratio in the computation of moments