  We present lattice QCD calculations of valence parton distribution function
(PDF) of pion employing next-to-next-leading-order (NNLO) perturbative QCD
matching. Our calculations are based on three gauge ensembles of 2+1 flavor
highly improved staggered quarks and Wilson--Clover valance quarks,
corresponding to pion mass $m_\pi=140$~MeV at a lattice spacing $a=0.076$~fm
and $m_\pi=300$~MeV at $a=0.04, 0.06$~fm. This enables us to present, for the
first time, continuum-extrapolated lattice QCD results for NNLO valence PDF of
the pion at the physical point. Applying leading-twist expansion for
renormalization group invariant (RGI) ratios of bi-local pion matrix elements
with NNLO Wilson coefficients we extract $2^{\mathrm{nd}}$, $4^{\mathrm{th}}$
and $6^{\mathrm{th}}$ Mellin moments of the PDF. We reconstruct the Bjorken-$x$
dependence of the NNLO PDF from real-space RGI ratios using a deep neural
network (DNN) as well as from momentum-space matrix elements renormalized using
a hybrid-scheme. All our results are in broad agreement with the results of
global fits to the experimental data carried out by the xFitter and JAM
collaborations.
