We present the first numerical investigation of the method proposed in Ref. [1] to utilize gradient flow to obtain precise determinations of higher moments of PDFs from lattice QCD, circumventing power divergent mixing with lower dimensional operators. We apply this method to obtain moments of the isovector PDF of the pion using four Stabilized Wilson Fermion ensembles with $m_{\pi}\simeq 411~\text{MeV}$ and lattice spacings $a\simeq 0.064, 0.077, 0.094$, and $0.12$ fm. We present preliminary results of ratios of three-point functions as a function of flow time, which can be used to extract the ratios $\langle x^2\rangle/\langle x \rangle$ and $\langle x^3 \rangle/ \langle x \rangle$.
We find that a significantly higher precision can be achieved with this method compared to the canonical approach, which requires boosting and cannot reach higher than the $\langle x^3 \rangle$ moment.