Pion distribution amplitude at the physical point using the leading-twist expansion of the quasi-distribution-amplitude matrix element

We present a lattice QCD determination of the distribution amplitude (DA) of the pion and the first few Mellin moments from an analysis of the quasi-DA matrix element within the leading-twist framework. We perform our study on a HISQ ensemble with $a=0.076$ fm lattice spacing with the Wilson-Clover valence quark mass tuned to the physical point. We analyze the ratios of pion quasi-DA matrix elements at short distances using the leading-twist Mellin operator product expansion (OPE) at the next-to-leading order and the conformal OPE at the leading-logarithmic order. We find a robust result for the first non-vanishing Mellin moment $\langle x^2 \rangle = 0.287(6)(6)$ at a factorization scale $\mu=2$ GeV. We also present different Ans\"atze-based reconstructions of the $x$-dependent DA, from which we determine the perturbative leading-twist expectations for the pion electromagnetic and gravitational form-factors at large momentum transfers.


I. INTRODUCTION
The study of inclusive deep-inelastic processes by describing them using process-independent parton distribution functions (PDFs) has resulted in a good understanding of collinear internal structures of hadrons. The next generation of experimental facilities (e.g., [1][2][3][4]) will focus on observables that characterize hard semi-inclusive and exclusive processes to relate intrinsic properties of hadrons to those of the partons. The generalized parton distribution functions (GPDs) and the meson distribution amplitudes (DAs) will play crucial roles as universal soft-functions in the descriptions of exclusive reactions in hard kinematical regimes through factorization. Thus, non-perturbative determination of such parton distributions and amplitudes are currently essential.
Concretely, the pion DA, φ(x), captures the overlap of the pion with a state with two collinear valence quarks carrying fractions x and (1 − x) of the pion light-front momentum P + [5][6][7]. Hence, the pion DA is of theoretical interest due to its proximity to being the lightfront wave-function [6] of the Nambu-Goldstone boson of chiral symmetry breaking; by comparison with DAs of non-Goldstone pseudoscalar mesons, one could learn about how the fundamental quark-gluon interaction leads to the special properties of the pion (see review [8], and Ref. [9] for our lattice calculation in a related direction). Phenomenologically, the properties of the pion DA, such as its shape and its Mellin moments, are still not precisely known and the main experimental input has been * gaox@anl.gov † nkarthik.work@gmail.com ‡ yong.zhao@anl.gov from the factorization of the pion-photon transition form factor [10][11][12][13][14][15], and from past [16][17][18][19][20][21][22][23] (and ongoing [24]) investigations of the large-Q 2 behavior of the pion electromagnetic form factor [5,7,25]. From a field-theoretic standpoint, the determination of the pion DA involves the light-front correlation [5,7,26], with a dimensionless invariant amplitude I renormalized in the MS scheme at scale µ defined as, if π P + I(λ, µ) = 0 d(−z − /2)γ + γ 5 W + u(z − /2) π + ; P , (2) where W + is a straight Wilson line from −z − /2 to z − /2 along the light-cone. Various model-based determinations (e.g., [27][28][29][30][31]) of φ have given key insights into the full x-dependence and the moments. For a rigorous QCD-based non-perturbative calculation, one needs to rely on lattice QCD computation. However, the unequal time-separation in the above light-front correlator has prevented a direct computation of DA.
The first few Mellin and Gegenbauer moments of the pion DA have been computed from leading-twist local operators on the lattice [32][33][34][35][36][37][38][39][40][41][42]. The difficulty of this approach is the nontrivial renormalization of these operators on the lattice, which limits the number of lowest calculable moments. One way to avoid this challenge is through a leading-twist expansion of the pionto-vacuum transition matrix element of certain spatially separated operators [43]. The short-distance logarithmic divergences are absorbed as part of perturbatively computed Wilson coefficients, and most importantly, the expansion coefficients are proportional to the moments of the DA. Such a leading-twist expansion approach was first applied using the pion-to-vacuum transition matrix element using two current operator insertions [43], which has been been utilized in further studies of DA and the moments [44,45]. Analogous real-space analysis using current-current correlators was also applied to the pionto-pion forward matrix element to determine the pion PDF [46][47][48]. Besides, there are also approaches based on the leading-twist expansion of current-current correlators in the Fourier space [49], including the method that uses an intermediate heavy quark [50][51][52], to calculate the higher moments of DAs and PDFs.
Recently, there have been new advancements in the determination of parton distribution functions using multiplicatively renormalized, spatially-extended quarkantiquark operators based on the LaMET [53][54][55] and the pseudodistribution [56][57][58] approaches. For example, Refs. [59][60][61][62] are recent computations of the pion PDF using the bilocal quark-antiquark operator. For DA, one can construct the pion-to-vacuum transition matrix element of such a bilocal quark-antiquark operator [63], which we simply refer to as the quasi-DA matrix element. The aim of this paper is to apply the leading-twist expansion of the quasi-DA matrix element in a manner similar to Ref [43], to extract the Mellin moments of the pion DA, and attempt to reconstruct the shape of the pion DA based on strategies utilized in the case of PDFs. Previously, such quasi-DA matrix elements of both the pion and kaon have been investigated using the x-space LaMET matching [64][65][66][67][68]. Due to the differences in the renormalization method and the analysis methodology for a leading-twist expansion approach, we expect this work to shed new light into the quasi-DA method as a way to obtain the pion DA and its moments. Apart from the above leading-twist expansion or effective theory matching approaches, it has also been proposed to directly obtain the structure functions from the hadronic tensor on the lattice through a nontrivial inversion of Euclidean correlators [69], which has not yet been applied to the extraction of DAs.
The plan of the paper is as follows. First, we give an overall description of our methodology in Section II. Then, in Section III, we present the perturbative results pertaining to conformal and Mellin OPEs. In Section IV, we give the details of our lattice calculation. In Section V, we present the details of the determinations of the bare quasi-DA matrix element and the renormalized ratios thereof. In Section VI, we present the results on the pion DA; here, we first present the analysis specifications, then we present the Mellin moments for fixed spatial distances, after which we present our model-independent determination of the first two Mellin moments, and finally, we describe the model-dependent reconstructions of the shape of the pion DA. From such reconstructed xdependent DA, we present the perturbative expectations for pion form-factors at large momentum transfers. In Section VII, we summarize our findings.

II. METHOD
We specify four-vectors as v ρ , whose components are (v 0 , v 1 , v 2 , v 3 ) with v 0 being the temporal component, and v = (v 1 , v 2 , v 3 ) as the spatial component. We use the ρ = 3 direction for spatial separations and as the direction of the pion momentum. The metric convention is v · w = v 0 w 0 − v · w. We specify the Dirac γmatrices as γ ρ , and we use the Minkowskian convention for them. For the ease of understanding, we specify them as (γ t , γ x , γ y , γ z ) respectively, when explicitly mentioning a matrix. We specify the bare and renormalized quantities using superscript "B" and "R" respectively.
We use the short-distance behavior of the quasi-DA matrix element of a boosted pion, π + (ud), to determine its leading-twist DA. The quasi-DA operator is the equaltime bilocal quark bilinear operator, with the straight Wilson-line W −z/2,z/2 connecting the quark and antiquark that are separated spatially as z = (0, 0, 0, z 3 ). At non-zero z, the operator suffers from a linear divergence in the self-energy of the Wilson-line, e −c|z| , and also from the end-point logarithmic divergence, and therefore, it needs to be renormalized. Thus the operator above is bare, and hence, the superscript B to specify this. Let O R ρ (z, µ) be the renormalized operator in the MS scheme at scale µ. The quasi-DA matrix element for the pion is, with the on-shell pion momentum P = (E(P 3 ), 0, 0, P 3 ). The Lorentz invariant λ = −z · P = z 3 P 3 is called the Ioffe-time or light-cone distance in the literature. We note that the left-hand side of the above equation is not a Lorentz decomposition, instead we have defined h above in a form convenient for the leading-twist expansion that is proportional to P ρ 1 . As a specific case, z 3 = 0 and for all momentum P 3 , the local operator O R ρ (0) is the axialcurrent operator, and h(z 3 = 0, P 3 ) = f π , the pion decay constant.
The idea used in this paper is that for quark-antiquark separations z 3 that are small in QCD scales and for momenta P 3 > 0, one can describe the λ and z 2 dependencies of h R (λ, z 2 , µ) within a leading-twist OPE framework valid up to higher-twist contributions. This lets us relate the lattice-calculable equal-time quantity h R (λ, z 2 ) to the light-cone distribution amplitude φ(x, µ) at a factorization scale µ and its Mellin moments, 1 It must be implicitly understood that h is also labeled bŷ nρ = zρ/ √ −z 2 , as will be reflected in the ρ-dependent Wilson coefficients (c.f. [70]) in the operator product expansion (OPE) of O R ρ (z), and hence, of the leading-twist expansion of h R (z).
The framework is similar to the one used in the determination of the parton distribution functions (PDFs), however, the key difference for the case of DA is that the matrix element in Eq. (4) is between the boosted pion state and vacuum state. This results in a different leading-twist expansion than in the case of the forward matrix element for PDFs. As we will show in the paper, the leading-twist expansion of h R (λ, z 2 , µ) for DA is where C n,m (µ 2 z 2 ) are the Wilson coefficients calculable in perturbation theory that relates h to the DA via its Mellin moments at scale µ. By the superscript "tw2" we mean that the expansion ignores all terms with twists bigger than two. Henceforth, we will refer to Eq. (6) as the Mellin OPE (M-OPE). In Section III, we provide the NLO results for C n,m for O R 3 . Under an ERBL [5][6][7] evolution in µ, the different Mellin moments mix, which is also reflected in the non-vanishing off-diagonal nature of C n,m . At the level of leading logarithms and up to finite O(α s ) corrections, massless QCD is conformal and this helps in diagonalizing the leading-twist expansion (see review [71]) with respect to evolution. Such an expansion in terms of the conformal partial waves F n (λ/2, µ 2 z 2 ; α s ) is given as where a n are the Gegenbauer moments at scale µ, and they satisfy the simpler LO DGLAP evolution. In Section III, we provide the expressions for F n . We will refer to Eq. (7) as the conformal OPE (C-OPE). At next-to-leading order, the Mellin and Conformal OPEs differ in the finite α s terms and are the same up to α 0 s and α s log z 2 µ 2 terms. At the level of practical implementation, we have to truncate Eq. (6) and Eq. (7)for C-OPE, it is a truncation in conformal partial waves that each have infinite order in λ, whereas for M-OPE, it is a truncation in the order of λ. We will use M-OPE primarily in this paper to implement the matching at NLO, and compare the results with that obtained with C-OPE to cross-check that the results are approximately the same. The expressions in Eq. (6) and Eq. (7) are general for any pseudoscalar mesons, but it gets further simplified for the pion; due to isosopin symmetry, the Mellin and Gegenbauer moments for odd n vanish, and therefore, at leading-twist the matrix elements are purely real. We will use this fact in our analysis and set odd n moments to zero. Since the Gegenbauer moments for all n > 0 approach zero under evolution to µ → ∞, the C-OPE expression simply approaches F 0 (λ/2) asymptotically.
In the above discussion, we assumed that the operator is renormalized in the MS scheme, which cannot be directly implemented on the lattice. Furthermore, it was shown in Refs. [72,73], that the operator O B 3 (z), that has the γ z γ 5 structure, is multiplicatively renormalizable, whereas the choice O B 0 (z) that has the γ t γ 5 structure, mixes with theūγ 5 γ z γ t d operator when lattice-regulated fermions that break chiral-symmetry are used. Therefore, in this work, we only work with the O B 3 (z) operator to avoid the mixing. We adapt the renormalization group invariant (RGI) ratios [57,60,74] of hadronic matrix elements for the renormalization. Since the renormalization of O B 3 (z) is purely multiplicative, the ratio, with respect to the matrix element at a fixed momentum P 0 is an RGI quantity that we can determine on the lattice. We obtain the leading twist expression for M from the MS expressions for h tw2 in Eq. (6) and Eq. (7) by making use of its RGI nature as The actual lattice data in the range of z 3 and P 3 that we use could suffer from lattice corrections and higher-twist corrections to the continuum leading-twist expressions in Eq. (6) and Eq. (7). We model the two corrections using some functions L(z, P, a) and H(z, P ) respectively. With such corrections, we use expressions of the type, to perform our fits to the lattice QCD data for M(λ, z 2 , P 0 ) to obtain information on the Mellin or Gegenbauer moments of the DA depending on whether M-OPE or C-OPE is used for h tw2 , respectively. Equivalently, by modeling the functional form of φ(x, µ), we also reconstruct the x dependence of the pion DA. Alternatively, from such analyses, we can also infer the MS light-front Ioffe-time distribution (ITD) as We discuss the implementation of the above set of steps further in the following sections presenting our results.

III. ANALYTICAL RESULTS FOR CONFORMAL AND MELLIN OPE
A. Short distance factorization of the quasi-DA matrix element In x-space, the quasi-DA can be perturbatively matched onto the light-cone DA at large momentum, where the matching coefficient has been derived in the MS scheme at one-loop order [66]. In coordinate space, the quasi-DA matrix element h(λ, z 2 , µ 2 ) can also be perturbatively matched onto the light-cone correlation I(λ, µ) through a short-distance factorization formula, which has been derived in QCD at one-loop order [75] as where the matching kernel The 2δ ρ3 term in the curly bracket can be inferred from the factorization of the forward quasi-PDF matrix elements [70].

B. Conformal OPE
The LCDA can be expressed as the sum of Gegenbauer moments, where C 3 2 n (x) is a Gegenbauer polynomial (refer [76, Table 18.3.1]). The Gegenbauer moment a n (µ) can also be projected from φ(x, µ) as a n (µ) = 4(n + 3/2) 3(n + 1)(n + 2) At leading logarithmic (LL) accuracy, QCD is conformal, and φ n (µ) evolves multiplicatively with the anomalous dimension with H n = n i=1 1/i. The value of γ n is different from that for the Mellin moments of PDFs by a minus sign. Therefore, the Gegenbauer moments should be the basis of OPE under the conformal approximation to QCD.
The conformal OPE of the quasi DA matrix element in the MS scheme can be inferred from that for the currentcurrent correlator [43] as where λ = zP z , and the LL resummed coefficient with Γ and J n being the standard gamma function and Bessel function of first-kind respectively, and c n = 1 + α s C F 2π 5 + 2n 2 + 3n + n 2 + 2δ ρ3 2 + 3n + n 2 which is the same as the Wilson coefficients in the OPE of the helicity quasi PDF matrix elements [70], and which is the anomalous dimension of the nonlocal operator O ρ (z, µ). Note that in Eq. (18) the running of strong coupling is turned off because we assumed conformal symmetry.

C. OPE in terms of Mellin moments
If we do not include the scale evolution in the OPE, then we can consider expansion in terms of the Mellin moments in Eq. (6). The coefficient functions can be obtained by the relation where C (1) n,m can be read off from the coefficients of x m . The lowest few coefficient functions are L .

IV. COMPUTATIONAL SETUP
We used a mixed fermion action setup consisting of a 2+1 flavor HISQ sea quark action and a clover-improved Wilson fermion action for the valence quarks. The HISQ ensemble [77] was generated by the HotQCD collaboration, and consists of L 3 s × L t = 64 3 × 64 lattice sites at a lattice spacing of a = 0.076 fm. The sea quark mass in the setup corresponds to a near physical pion mass of 140 MeV. The tadpole improved Wilson clover valence quarks couple to 1-HYP smeared gauge links [78]. We tuned the Wilson mass to obtain the valence pion mass of 140 MeV. Therefore, both the sea and valence quarks are tuned to the physical point.
In order to compute the quasi-DA matrix element of a pion with momentum P 3 , the two essential ingredients are the π-π and π-O 3 correlators . The pion-pion twopoint function at source-sink time separation of t s is where with π † (P, t s ) = x π † (x, t s )e iP·x . The u s and d s represent Coulomb-gauge Gaussian smeared quark operators, with the smearing radius as 0.59 fm. At nonzero spatial momentum P = (0, 0, P 3 ), we implemented the boosted quark smearing [79] with quark boosts, k 3 = ±ζP 3 , for u and d respectively. The other ingredient, the pion-quasi DA-operator correlator with time separation t s , is The spatial part of the quark-antiquark separation where U HYP is the 1-HYP smeared gauge link, the same as those used in the Wilson Dirac operator. The u and d quark operators in Eq. (32) are not Gaussian smeared. One should note that the coordinates of the antiquark and quark are at x and x + z which differs from the one in Eq. (3), and hence the tilde on top of O 3 to make this distinction clear. This is due to the ease of implementation of the former convention on the lattice, and we defer the conversion to the analysiswise convenient convention in Eq. (3) at a later stage by multiplying the results with a phase exp (−iP 3 z 3 /2). We computed the quark propagators that occur in the Wick contractions of Eq. (29) and Eq. (31) on GPUs using the multigrid algorithm [80] as implemented in the QUDA suite [81][82][83].
We performed the above set of computations at eight different spatial momenta P = (0, 0, P 3 ), for n 3 ∈ [0, 7]. Thus, the highest momentum we use in this work is 1.78 GeV which is sufficiently larger than Λ QCD , and at the same time corresponds to P 3 a = 0.69 which is below the lattice-like scales. For these momenta, we decided to choose the phase parameter ζ in momentum smearing such that ζn 3 = 2 for n 3 ≤ 3 and ζn 3 = 5 for n 3 > 3, so that we could reuse smeared sources for multiple momenta to balance the computational cost and the signal-to-noise ratios.
We used 350 statistically independent configurations. We effectively increased the statistics many folds using the all-mode averaging method [84], implemented using exact inversions of the Dirac operator at N ex source locations x 0 , and sloppy inversions at N sl source locations. For n 3 ≤ 3, we used (N ex , N sl ) = (4,80), and for n 3 > 3, we used (8, 160).

V. DETERMINATION OF MATRIX ELEMENT
In this section, we discuss the details of the determination of the ground-state matrix elements of the bare quasi-DA operator at different momenta and quarkantiquark separations, and the RGI ratios that we construct from them.

A. Bare matrix element
We used the spectral decomposition of C ππ (t s ; P 3 ) and C πÕ3 (t s ; P 3 , z 3 ) to extract the bare quasi-DA matrix element. That is, the pion-pion correlator, and the pion-quasiDA correlator, where the kets are relativistically normalized, and Z n = π; P 3 π † (P 3 ) 0 , which we assume to be real and positive in this work. The summations in Eq. (34) and Eq. (35) run over all the eigenstates at definite momentum P 3 . However, for practical considerations, one truncates them including only the lowest N st states. We refer to the N st = 2 truncated expression as the two-state ansatz, and the N st = 3 truncation as the three-state ansatz. The periodicity of the two correlators is imposed above; for C πÕ3 (t s ) the periodicity can be seen by reflec- , along with q → γ z q,q →qγ z where q is either u or d, which is a symmetry of the Euclidean path integral, and therefore, C πÕ3 (t s ) = C πÕ3 (−t s ). We obtained the bare quasi-DA matrix element from the analysis of the ratio, whose spectral decomposition is simply the ratio of Eq. (35) and Eq. (34). It can be seen that the leading term for large t s behaves as R(t s ) → P 3h B (z 3 , P 3 )/Z 0 . First, we obtained the best fit values of the spectral parameters E n and Z n from the analysis of C ππ correlators. In a previous work [85], we discussed our fits to the pion correlator on the same gauge ensemble. In this work, we used the two-state and three-state fits with the ground-state energy fixed to the continuum dispersion, E 0 (P 3 ) = P 2 3 + m 2 π . We chose the fit ranges t s ∈ [t min , t max ] for C ππ (t s ) such that they covered the range used for the subsequent fits to the ratio R to be discussed next. Namely, we used t s ∈ [4a, 32a] for two-state fits and t s ∈ [2a, 32a] for three-state fits. In this way, we obtained good effective values of E n and Z n that best describe the excited state contribution to the ratio R(t s ) in the range of t s we made use of. However, as observed in Ref. [85], the values of E 0 , E 1 and E 2 from our final fit choice were within errors of the results when larger t min were used, albeit with noisier determinations. Whereas the value of E 1 at P 3 = 0 is consistent with the pole mass of π(1300), the value of E 2 at P 3 = 0 from three-state fits is much higher than expected at 3 GeV. Thus, as noted in [85], it is likely that the three-state fit with E 2 capturing the tower of excited states above E 1 via a single effective state.
In the next step, we used (Z n , E n ) from two-state and three-state fits on jackknife samples as inputs in our fits to R(t s ) over ranges t s ∈ [t min , t max ] on the same jackknife samples. The fits for R(t s ) used the above spectral decomposition with fit parameters being the amplitudes 0|O 3 |E n , and therefore, the fits were linear. We performed these fits to R(t s ) using two-state and threestate ansatz. For two-state fits to R, we chose t min = 6a, whereas for the three-state fits, we used t min = 4a. For both two-and three-state fits, we chose the maximum range of the fits t max = 20a for momenta n 3 ≤ 3 and t max = 15a for n 3 > 3 to avoid noisier estimates at larger t s . In this way, we extrapolated the ratio to t s → ∞ to obtain the bare quasi-DA matrix elementh B (z 3 , P 3 ).
In Fig. 1, we show some examples from our extrapolations of the real part of R. From top to bottom, the data are from momenta n 3 = 1, 3, 5 and 7 respectively. Let us first focus on the left panels. We show the data for Re(R) as a function of t s for few sample values of z 3 as specified near the data points. The magnitude of R is not important as it still depends on the two-point function amplitude Z 0 , and only its t s dependence is important here. For n 3 = 1 to 3, the variations with t s is smaller due to larger energy gap, E 1 −E 0 , which is about the gap between pion mass and that of π(1300). For larger n 3 , the variation of R with t s is significant due to the states being relativistic. Thus, the extrapolations using spectral decomposition of R is necessary in our calculation, especially in the important large momenta data set. The blue and the black bands show the extrapolations using the two-state and three-state ansatz respectively. Within the statistical errors, the two extrapolation bands satisfactorily describe the t s dependence of the lattice data. However, the three-state fits have a tendency to be closer to the central values of the data when compared to the two-state ones. In the right panels of Fig. 1, we compare the resultant values of the bare quasi-DA matrix element,h B (z 3 , P 3 ), from the two-state and three-state The pion decay constant fπ, modulo the finite renormalization factor ZA, is shown as a function of momenta P3 = 0.254n3 that is used in the extraction. The results using Γ = γzγ5 and γtγ5 are shown in the plot. The dashed curve is the value of fπ/ZA from n3 = 1.
fits to R as a function of z 3 . Note that the ordering of blue and black points in the right panel is not in one-toone correspondence with the left panel as there is also an additional factor Z 0 that is different between the left and right panels. Within errors, the two-state and threestate extrapolated values are consistent, with perhaps a slight tension in the n 3 = 3 momenta. The errors on the three-state fits are comparable or smaller than in the two-state fits due to the t min being 4a for three-state ones compared to 6a for two-state fits. The relative error at larger momenta n 3 > 5 increases at even shorter z 3 . Nevertheless, as we will discuss in the next subsection, the growth in statistical error in shorter z 3 is reduced due to the RGI ratios that we will construct, and, due to the same reason, even the slightest discrepancies between the two-state and three-state fits seen in the bare matrix element in Fig. 1 will be reduced further. We found similar consistency between the two-state and three-state fits in the case of Im(R). The z 3 = 0 value ofh B is the bare pion decay constant, f π /Z A , where Z A is the finite renormalization constant for the axial current operator. Thush B (0, P 3 ) has to be constant with P 3 if there were no systematical errors in the extrapolations and if O(aP 3 ) lattice corrections did not affect the lattice results, and therefore, provides a cross-check on our calculation. In Fig. 2, we show f π /Z A as a function of momentum P 3 used in h B (z 3 = 0, P 3 ). The red circular data points are the results using the O 3 (z = 0) operator. In addition, we also looked at the local matrix element from O 0 (z = 0). We show those values as the blue triangles in Fig. 2. The values of f π /Z A are consistent with being constant with respect to P 3 , with perhaps a slight dip in the central value around n 3 = 3, which could be due to statistical fluctuation. The excited state contribution and the Lorentz structure of the O 3 and O 0 matrix elements are different, and hence, the consistency between the f π /Z A determinations from the two observables is reassuring. Using RI-MOM renormalization procedure, we determined Z A = 0.969(1) for the ensemble used in this paper (see Appendix A). From the most precise values of the matrix element at n 3 = 0 for O 0 (z = 0) and n 3 = 1 for O 3 (z = 0), the values of f π from the two observables are 130.0(4) MeV and 129.7(4) MeV respectively. These results agree with the FLAG average for 2+1 flavor QCD f π = 130.2(8) MeV [86]. Due to the normalization condition x 0 = 1 that we will impose on the matrix elements, f π will not play any further role in this calculation.

B. Renormalized ratios
We used the RGI ratios of h B (z 3 , P 3 ) to get the renormalized quantities. First, we shifted the location of the operatorÕ(z) by −z/2 in order to conform with the definition in Eq. (3). We did this by multiplyingh B (z 3 , P 3 ) with a phase exp(−iz 3 P 3 /2) from the translation. Next, we improved the ratio in Eq. (8) to impose the condition that the ratio should be exactly unity at z 3 = 0 using the so called double ratio procedure. Thus, in the end, we determined the RGI ratio [57,60,74] as We have written the arguments in the right-hand side above in terms of z 3 and P 3 to make the definition clear. The factor in the second parenthesis above should be exactly one, devoid of any systematical and statistical errors. We refer to the fixed momentum P 0 used to form the ratio as the reference momentum. We used a nonzero value of P 0 for two reasons -1) the leading-twist part of the matrix element of O 3 vanishes at zero spatial momentum. 2) by using larger P 0 , the higher-twist corrections, such as (Λ 2 QCD z 2 3 ) k , present in h R (z 3 , P 0 3 ) are made relatively smaller compared to the leading-twist terms containing powers of P 0 3 z 3 . However, if we use P 3 < P 0 3 , the possible advantage of larger P 0 3 is rendered meaningless. Therefore, we only used P 3 > P 0 3 . In Fig. 3, we show the real and imaginary parts of the RGI ratio with reference momentum P 0 3 = 0.254 GeV. In the left panel, we show Re M as a function of z 3 for easier visibility of data at different P 3 . We show the results for Re M obtained using the bare matrix elements from the two-state and three-state fits together. It is clear that the two extrapolated results are quite consistent with each other, even more so after forming the RGI ratios, wherein any correlated systematical errors could get canceled between the numerator and denominator. It is also striking that the errors at larger momenta at small to moderate range of z 3 are statistically well determined, thanks to the statistical correlation in the data at P 3 and P 0 3 at a given z 3 . By construction, the z 3 = 0 value of M is exactly 1. In the rest of the paper, we will use the data for M obtained using three-state fits. In the right panel of Fig. 3, we show Im M from three different representative values of P 3 . Due to chiral symmetry, the leading-twist part of h R (λ, z 2 , µ), and hence M, should be purely real. This is demonstrated in our data by the vanishing of Im M well within statistical errors at different P 3 and z 3 . Therefore, we only analyzed ReM and imposed the symmetry of pion DA about x = 0 explicitly by setting x n = 0 for odd n.

A. Analysis strategy
We extracted the leading-twist DA related quantities from Re M by fits to the corrected (as well as the uncorrected) leading-twist expression Re M tw2,corr in Eq. (10).
Let P be the set of free parameters that enter Re M tw2,corr ; for example, they could be the set of moments or the parameters of a DA ansatz. We found the best fit values of P by the standard χ 2 fits using ]. We chose z min 3 > a to reduce the effect of lattice corrections at lattice-like separations. We used z max 3 =0.456 fm, 0.608 fm and 0.76 fm to take into account possible variations in the fitted values due to higher-twist contaminations that we did not capture in Eq. (10), and at the same time remain in moderately small values of z 3 that are allowed given the constraint of the lattice spacing we are using. We used the momenta from P min 3 > P 0 3 , the reference momentum. We used the full covariance matrix Σ to take care of correlations between the data at different z 3 and P 3 . The value of the strong-coupling constant α s enters the leading-twist OPE expressions. At NLO, the scale at which it needs to be determined is ambiguous. In this work, we use the value of α s (µ) at the same scale at which the DA is determined, namely, at µ = 2 GeV. We take the value of α s (2GeV) = 0.303 determined from the running of α s taken from the PDG [87].
The leading-twist expansion approach comes with systematic uncertainties due to the possible analysis choices, such as, the choices of z min/max 3 , P min/max 3 , and the choice of ansatz for higher-twist corrections, to list a few. Apart from presenting a scatter of the fitted results for all possible combination of analysis choices, it is helpful to summarize a result compactly to capture its central value, the systematic spread due to analysis variations, and the statistical error on the central value. To achieve this, we used the following procedure. Let S be the set of analysis choices, and let P a be the set of best fit parameters for a particular choice a ∈ S. Following the approach presented in [60,88] closely, for some function of parameters, F (P), we first found the mean valueF and standard deviation σ width of results for F over all analysis choices in a given jackknife block as for weights w a for each analysis choice. From the central value and width of scatter per jackknife sample, we determined the final estimate of F as (mean)±(statistical error)±(systematic error) by finding J av (F ) ± J er (F ) ± J av (σ width ), where J av (. . .) is the jackknife average and J er (. . .) is the jackknife standard error of a quantity. One possibility for the w a is the Akaike information criterion (AIC) weight given by e − 1 2 (χ 2 +2f ) where f is the number fit parameters. We found that such an estimator for our case had a tendency to choose only a few of the analysis choices and does not represent the true scatter present in our analysis. Therefore, we followed the approach we used in our earlier work [60], which is to set a constant w a (i.e., unweighted averaging) for all analysis choices. In this way, we took all the choices on equal footing.
B. Fixed-z 2 analysis: looking for corrections to continuum leading-twist expectation The simplest analysis of the data is to study the λ dependence of M(λ, z 2 3 ) at fixed values of z 3 (refer [89] where the idea was first proposed for the forward matrix element). In this way, the variation in λ comes only from the variation in P 3 . The output of this analysis is the set of Mellin or Gegenbauer moments at a fixed scale µ, based on whether M-OPE or C-OPE is used respectively, as a function of z 3 . Using the degree of agreement of the z 3 -dependent moments with a plateau in z 3 is a nice way to understand whether the fixed-order leading-twist framework is applicable to the lattice data in a range of z 3 or not, and which type of corrections to continuum leading-twist expansion are seen. Such a diagnosis of the lattice data has been performed in the case of the forward matrix element in the PDF determination [60]. Here, we apply such an analysis for the off-forward matrix element for the first time.
We used the purely leading twist expression (i.e., we set the higher-twist correction H and lattice correction L to zero) in Eq. (9). By using the C-OPE for the leadingtwist expression, we obtained the Gegenbauer moments a n by using them as the fit parameters. Similarly, by using M-OPE, we obtained the Mellin moments x n . Since each Mellin or Gegenbauer moment adds an additional fit parameter, we needed to truncate the OPE at a finite order 2N max which contains N max number of even-n moments; we successively increased N max from 2 to 4. For C-OPE, the fits became unstable for N max > 3 as the dependence on Gegenbauer moments beyond a 2 is rather weak. With M-OPE, we were able to go up to N max = 4 in this analysis.
In the top panel of Fig. 4, we show x 2 as a function of z 3 upto z 3 = 0.91 fm by applying M-OPE to M with P 0 3 = 0.254 GeV. We show the results using three different truncations N max . For z 3 < 0.8 fm, we can see that N max = 4 is sufficient given the data precision. The near plateau in x 2 around a value of about 0.28 shows that the leading-twist expansion describes the lattice data for M to a good approximation and that any higher-twist or lattice corrections are subdominant. This provides a validation of the leading-twist fixed-order perturbative framework that we are using to describe the nonpertubative lattice QCD data in the range of subfermi values of z 3 we use. The near plateau in x 2 around a value of about 0.28 shows that our fitting form based on leadingtwist perturbative expansion at NLO, with the choice of µ = 2 GeV, can describe the lattice data within the current statistical error. However, as shown in Refs. [62,90], perturbation theory may become unreliable at large z due to the resummation of large ln z 2 µ 2 [91], so a more dedicated study of the comparison between fixed-order and renormalization-group improved OPEs needs to be done to understand the results we have observed. Nevertheless, the central value of x 2 changes by about 0.01 (≈ 4%) as z 3 is increased from 0.1 fm to 0.6 fm. Hence, in the subsequent analysis of the data, we will include correction terms such as H and L to the leading-twist expansion. Since the correction itself is small, we were not able to deduce a possible functional form for the functions H and L that might be present, as we did in the case of the pion PDF in Ref [60]. In the middle panel of Fig. 4, we show a similar z 3 dependence of x 4 at µ = 2 GeV. Since the analysis only made use of six different data points at each z 3 , the errors on x 4 are larger. Significant information on x 4 enters M(λ, z 2 ) only beyond z 3 > 0.5 fm, wherein we find initial indication that x 4 ≈ 0.15, as we will find in the subsequent combined analysis of all the data. Within the larger errors, the data is consistent with z 3 independence.
In the bottom panel, we compare the values of x 2 obtained using M-OPE (black squares) with that from C-OPE (red circles). For C-OPE, we converted the values of the fitted Gegenbauer moment a 2 to x 2 through the simple linear relation (e.g., [92]), x 2 = 1/5 + 12/35a 2 . While results from both M-OPE and C-OPE are approximately z 3 independent, the values of x 2 from M-OPE is about 3% higher than that from C-OPE. Since both M-OPE and C-OPE have converged well with respect to the OPE truncation, this is likely due to the remnant finite O(α s ) corrections that are missing from C-OPE, but captured correctly in M-OPE. We also show the result using α s = 0 in the M-OPE, which we refer to as the tree-level result. Surprisingly, the tree-level result is also approximately plateaued, showing the effect of perturbative ln µ 2 z 2 3 in the OPE to be mild in the range of z 3 that we investigated using the choice of scale µ = 2 GeV. Henceforth, we will primarily show the results using M-OPE, and use C-OPE to compare those results with and serve as an indirect way to quantify the perturbative uncertainty by going from leading-log order to NLO.

C. Determination of Mellin moments
Having shown the effectiveness of leading-twist OPE in capturing the λ and z 2 3 dependencies of M(λ, z 2 3 ), we now comfortably apply the framework to estimate the lowest two Mellin moments x 2 and x 4 through a combined fit to all the data for M(λ, z 2 3 ) within a range of z 3 . The fits are similar to the ones in the last section with the moments themselves as the free fit parameters. Therefore, the method is independent of any modeling of the x-dependence of the DA.
In the absence of any obvious visible evidence in Fig. 4 for lattice and higher-twist corrections, we simply modeled them. For the lattice correction that affects the lattice-like separations z 3 , we assumed a possible presence of (P 3 a) 2 corrections in the quasi-DA matrix element similar to the one in the quasi-PDF matrix element [60]. Thus, we chose a functional form, with l 2 being a real valued fit parameter in the modeled correction. In this way, such a term can effectively affect the leading twist terms at O(λ 2 ) through a z −2 3 type correction to the moments. In addition, there could be momentum independent O(a) or O(a 2 ) corrections; since the continuum leading-twist expressions work quite well in describing the lattice data at a fixed lattice spacing, it is likely that such momentum independent corrections can be absorbed as part of the Mellin moments and the extracted DA themselves at that finite lattice spacing. For the higher-twist corrections, we followed a procedure similar to the one in [43], and assumed that the corrections resemble the one from twist-4 DA terms captured via a conformal OPE. For this, we added terms of the form, with h m being the free parameters. In the above equation, we used the tree level conformal partial waves F (0) m , obtained by setting α s = 0 in Eq. (18), to avoid modeling the logarithmic z 2 3 dependencies using extra fit parameters. Since we introduce the corrections as a ratio via Eq. (10), the term H can start at O((λ) 0 ) as the condition that M tw2,corr (λ = λ 0 , z 3 ) = 1 is satisfied by construction. Thus, we included F 0 ∼ (λ) 0 as the leading term to Eq. (40). We used N HT = 0, 1, 2 in order to keep the number of correction terms required to be small and at the same time take into account the sensitivity of the extracted results on the modeled higher-twist effects. However, one should note that there exists more complex approaches to model the functional form of H (e.g., see [93] that uses renormalon model) than the simpler functional parametrization that we adopt in this work.
In Fig. 5, we show the ratio M(λ, z 2 3 , P 0 z ) as a function of λ = P 3 z 3 , using only the lattice data with z 3 < 0.91 fm. The top and bottom panels are obtained using the reference momentum P 0 3 = 0.254 and 0.508 GeV respectively. We show the lattice data points from different P 3 > P 0 3 together in the two plots, and we differentiate between them by the colors and symbols used. The insets in Fig. 5 simply magnify the range λ < 2. We used the NLO Mellin OPE for the twist-2 contribution in Eq. (10), with and without the H and L correction terms that we discussed above to fit the lattice data for M. Since the usage of non-zero P 0 3 is not common in  ]. The red circles are obtained using Mellin moments as fit parameters without constraints, whereas the green circles are obtained using a positivity constraint on the pion DA. The variability comes from the number of higher-twist correction terms NHT, the number of lattice correction terms NLC, the reference momenta n 0 3 used in the ratio, and the ft ranges. The complete specification (NHT, NLC, n 0 3 , z min the literature, we note that at a given fixed λ, the z 3 dependence comes from two sources even at leading-twist; namely, the ln(z 3 ) dependence due to perturbative evolution and a polynomial dependence due to terms such as P 0 3 z 3 in the denominator of the ratio. This is the reason that at P 0 3 = 0.508 GeV, one finds a somewhat larger z 3 dependence than one would expect only from the perturbative logarithm. To be clear, the presence of additional P 0 3 z 3 type polynomial dependence on z 3 is not a disadvantage as such terms are captured within a leading-twist framework without any modelling, and the consequent spoiling of near universality with respect to λ = z 3 P 3 is not a practical issue from the point of view of fits. We truncated the OPE using N max = 4 number of even-n moments. The dashed curves in Fig. 5 are the central values of the best fit curves when N HT = 1 and L(z 3 ) terms are used as corrections, whereas the solid curves are obtained without any correction terms. In the example fit shown, we used a range z 3 ∈ [2a, 0.608 fm]. It is clear that the two types of fits work quite well in describing the lattice data. We found the fits with higher-twist correction terms to perform marginally better in terms of χ 2 /df, and this shows up in the tendency for the dashed curves in Fig. 5 to pass closer to the central values of the lattice data points. Taking the case with N HT = 1 and P 0 3 = 0.25 GeV shown in the top panel as a sample case to discuss the typical values of the fit parameters that we obtained in our fits, we find We see that our data sufficiently constrains only the lowest two even Mellin moments. The lattice correction term l 2 does not impact the fits, whereas the higher-twist correction term h 0 cannot be neglected. The value of h 0 = (81 MeV) 2 is in the ball-park of the value of highertwist correction we empirically found in the quasi-PDF matrix element in Ref [60]. However, its value is small compared to the expectation for the twist-4 0-th Gegenbauer moment based on QCD sum rules [43,94], namely, δ 2 π ≈ (300 MeV) 2 . Apart from the one case shown in Fig. 5, we also repeated the fits over the following 72 combinations of analysis choices: a) N HT = 0, 1, 2, (b) with and without a lattice correction term, N LC = 0, 1, (c) z min 3 = 2a, 3a to take short-distance lattice artifacts into account, (d) z max 3 = 0.456, 0.608, 0.76 fm to take variations from higher-twist effects into account, (e) P 0 z = 0.254, 0.508 GeV for variation from reference momentum. In Fig. 6, we have shown the determination of x 2 and x 4 from each of the choices (N HT , N LC , n 0 3 , z min 3 /a, z max 3 /a) as a data point (red circles). We specify the analysis choice to the side of each point. We found χ 2 /df < 1.6 for all the fits, and hence were acceptable. To make the trend in the fit parameters with respect to the added higher-twist terms visible, we have grouped the points in Fig. 6 in three sets with N HT = 0, 1, 2 as indicated by the horizontal dashed lines. From the systematic shift in the determined values of moments, there appears to be a small but non-negligible effect of adding a z 2 3 F (0) 0 (λ/2) term to the leading-twist OPE. The addition of one more term z 2 3 F (0) 1 (λ/2) seems to only make the fits noisier. Thus, given the precision of the data, usage of N HT = 1 seems to be sufficient. We show the scatter of the other fit parameters as well as the values of χ 2 /df in the various fits in Fig. 13 in Appendix B.
Using the analysis method we discussed earlier, we summarize the content of the red points in Fig. 5 as the following unweighted averages with the statistical and systematic errors: at a scale µ = 2 GeV. These summary estimates are shown in the vertical bands of Fig. 5; the inner band includes the statistical error only, whereas the outer one includes both the statistical and systematic error. It can be seen that the outer band covers most of the scatter due to the various choices, and hence is representative of our data. If we assume that the pion DA is positive at all x at µ = 2 GeV, then we can improve the stability of fits by imposing inequalities on the Mellin moments that follows from φ(x, µ) > 0, and subsequent derivatives of x n with respect to n at values infinitesimally closer to integer values. Namely, as we explain in [60], we obtain the inequalities (a) x n > x n+2 and (b) x n+2 + x n−2 > 2 x n . In the analysis, we imposed the two constraints through a change of variables x 2n ≡

Nmax i=n
Nmax j=i e −λj . We have shown the results of such constrained fits as the green circles in the two panels in Fig. 5. We find the resulting values of the two lowest moments to be well determined, especially as the number of fit parameters is increased when we set N HT = 2. In the case of x 4 , such a procedure results in more precise estimates compared to the unconstrained estimates shown as red circles. We find from this constrained analysis that (11) (20). (43) at µ = 2 GeV. Using the model-independent estimates of the Mellin moments themselves, we can reach a few conclusions. The values of these two Mellin moments in the large Q 2 limit of DA, φ(x) = 4(1 − x 2 )/3, are x 2 = 0.2 and x 4 = 0.0857. These differ quite significantly from the values we determined, and hence we can conclude that the pion DA at the physical point and at a scale of µ = 2 GeV differs from the asymptotic DA. In fact, since x 2 > 0.2, we can expect the x-dependent DA to be flatter compared to the asymptotic DA. The DA cannot be a completely flat DA [95], φ(x) = 1/2 which is characterized by x 2 = 1/3 ≈ 0.33, and x 4 = 0.2. We can consider another extreme case of a double humped Chernyak-Zhitnitsky (CZ) DA [27], φ(x) = 15(1 − x 2 )x 2 /4, at µ = 2 GeV, with Mellin moments as x 2 = 0.4285 and x 4 = 0.2381. These values are not compatible with the values we find. Instead, if we assume a simple oneparameter ansatz, φ(x) = N (1 − x 2 ) α , and solve for α using the value of x 2 = 0.2866, we find the exponent should be around α = 0.244. In the next subsection, we perform more elaborate fits to such Ansätze.
We performed a similar set of analyses using the C-OPE for the leading-twist contribution in Eq. (10). However, we were not able to obtain stable fits with a more complex N HT = 2 correction term to C-OPE without imposing any constraints on the Gegenbauer moments, which resulted in a spurious negative-valued a 2 at the cost of a large-valued higher-twist coefficient h 1 . Since it is likely due to overfitting of the data, we excluded this analysis choice. To summarize, using all other combinations of analysis choices, we found a 2 = 0.227 (18)(23) and a 4 = −0.16 (13) (30). These values correspond to Mellin moments, x 2 = 0.2779(63)(79) and x 4 = 0.121(15)(28), using the linear relations between a n and x n (e.g., [92]). It is reassuring that two ways of incorporating the leading-twist OPE result in similar values of the first two Mellin moments. It is also clear that when written using C-OPE, the essential non-vanishing contribution mainly comes from the a 2 Gegenbauer moment. The non-vanishing value of x 4 we find using the Mellin OPE, while being non-trivial information from the perspective of M-OPE, becomes trivial when expressed in terms of non-vanishing a 2 and a vanishing a 4 from the C-OPE perspective.
As a way of estimating perturbative uncertainty, we performed the above set of analyses at scales of µ = 4 GeV and √ 2 GeV, using α s (4 GeV) = 0.227 and α s ( √ 2 GeV) = 0.3607. We then perturbatively ran the estimated Mellin moments to the fixed scale µ = 2 GeV using the NLO implementation that incorporates the mixing among the Mellin moments (as discussed in Ref. [29]). Through this procedure, we found x 2 , x 4 at µ = 2 GeV to be [0.2822 (85)(47) Finally we compare our findings for the Mellin moments with some recent lattice QCD calculations at µ = 2 GeV. The work [39] using a dynamical QCD simulation and using the local twist-2 operator approach obtain x 2 = 0.28(1) (2). Another series of works from RQCD that culminated in Ref [42] using the local operator approach obtain x 2 = 0.240 (6)(2)(3)(2)  pansion method using current-current correlators [45] results in a scatter of values around x 2 ≈ 0.3. Using the quasi-DA matrix element as used in this work, but using LaMET x-space matching, the work [67] estimates x 2 = 0.244 (30) (20), and the most recent work [68] using the hybrid-renormalization method [90] estimates x 2 = 0.300 (41). Our result lies in the ball park value of previous estimates, but it is about 2.4-σ (including statistical and systematic errors in both works naively as the net error) larger from the estimate using the local operator approach in Ref [42]. In the future, we need to investigate the remaining systematical uncertainties in our work that we did not quantify, such as the effect of finite lattice spacing, and see if the tension between the values of Mellin moments obtained with two completely different methods, reduces or persists.

D. Prior-sensitive reconstruction of the pion DA
Our determination of the lowest two Mellin moments in a model-independent manner is the important result in this paper. However, within the framework of fitting phenomenology motivated Ansätze to the lattice data, we can reconstruct the x-dependence of the DA, φ(x). For convenience, we define the variable u via so that the DA has support from 0 to 1. In principle, once we know all the Gegenbauer moments from fits to C-OPE, or inferred from M-OPE, we can perform a modelindependent reconstruction using The caveat that all the moments a 2n need to be known makes such an approach not usable in practice; as we saw, the real-space quantity M(λ, z 2 ) converges in the accessible range of λ < 6 rapidly with respect to the number of Gegenbauer moments a n (as the main content of M can be summarized approximately with a value of a 2 ), whereas the corresponding convergence in u (or x)space is rather slow. The problem is easy to understand by considering a behavior φ(u) = N u α (1 − u) α . In the last section, from the value of x 2 , we expected α ≈ 0.25, which differs significantly from the leading term with α = 1 in Eq. (45). We can improve the convergence by using a complete basis that is orthonormal with respect to a weight function, w(u) = u α (1 − u) α , rather than the weight function w(u) = u(1 − u) that the Gegenbauer polynomials C 3/2 n are orthonormal with respect to. Such an idea was pursued in [30] using the Gegenbauer polynomial basis C α+1/2 n (1 − 2u), which we follow in this paper. To impose the evenness of φ(u) around u = 1/2, we restrict the functions to even n. That is, we expand, with s 0 = 1. The value of α describing the family of complete functions is arbitrary, but a usage of α that is close enough to the large/small-x exponent leads to a better convergence with respect to the truncation order N G . In Fig. 7, we show a specific example of the better convergence of an example DA, φ(u) = 1.47u 0.2 (1 − u) 0.2 , when expanded in a nearby C 0.9 n polynomial basis as compared to an expansion in C 3/2 n polynomials. We note that the polynomials C α+1/2 n (1 − 2u) are proportional to another complete basis, the Jacobi polynomials, P α,β n (1 − 2u) for α = β, that have been proposed [96] as a good choice in the analysis of PDFs even when α = β.
First, we determined the best fit values of the exponent α of the one-parameter Ansatz, that best describes the lattice data via Eq. (10) for each analysis choice that we described earlier. Essentially, the parameter α enters the fits through the αdependent Mellin moments. Therefore, unlike the modelindependent analysis of moments that we presented in the previous subsection, all the moments are now related through a single unknown parameter α. We truncated the Mellin OPE at order N max = 6 as before. For different analysis choices, we found the values of the α to lie in a range between 0.2 to 0.32.
In the next step, we generalized the parametrization for the DA by an expansion in C 1/2+α n as given in Eq. (46). We followed the approach in Ref [88] to slowly generalize from the one-parameter Ansatz above to more flexible ones using a complete basis of functions to capture the corrections to Eq. (47). We used the best fit  values of α from the one-parameter fits from the previous step to choose the basis, C 1/2+α n . Even though the values of α did not change much based on the analysis choices, we took care of using the corresponding α values for a given analysis choice. The fit parameters are the coefficients s n in Eq. (46), which enter via the Mellin moments or Gegenbauer moments (the implementation of fits can be made computationally faster by pre-evaluating the moments of Gegenbauer Polynomials 1 0 (2u − 1) n u α (1 − u) α C 1/2+α n du, from which the Mellin moments are obtained as linear combinations.) We imposed the external constraint on the allowed amount of fluctuations about φ 1−param through constraints on the expansion coefficients that |s n | δ. We realized this via Gaussian priors on s n added to the χ 2 using the central values and widths of the priors of all s n being 0 and δ respectively. One should however note that there is a priori no expected value for δ simply from general considerations. From practical considerations, we will present the reconstructions by imposing successively weaker constraints from δ = 0.05 to δ = 0.2. For even larger values of δ, we found the reconstruction to be very noisy and oscillatory.
In the panels of Fig. 8, we show the reconstructions of DA at µ = 2 GeV using prior widths δ = 0.05, 0.1 and 0.2 from top to bottom respectively. For the cases we show, we performed the fits over a range z 3 ∈ [2a, 0.61 fm] and using a reference momentum P 0 3 = 0.254 GeV. In each panel, we show the reconstructions based on M-OPE without any correction terms, with N HT = 1, and with (N HT , N LC ) = (1, 1). The changes due to such variations are small, especially the effect of N LC being negligible. From δ = 0.05 to 0.1, the effect of relaxing the prior of s n is primarily to increase the statistical error on the bands while closely agreeing with the one-parameter reconstruction. The reconstructed DA starts becoming slightly oscillatory and with larger error band when δ is relaxed to 0.2. Unlike the case of PDFs, which typically show a subdued prior and model dependence, we found the reconstructed DA to be sensitive to the prior that is applied. This is not surprising given that the essential content in our quasi-DA matrix element in the range of λ we used is a 2 , and the problem posed by such limited information in the DA reconstruction is well known in the literature. However, the use of the C 1/2+α n basis was useful to quantitatively and systematically reconstruct the pion DA that depends on the extent to which one allows the DA to deviate from the default one-parameter model. Thus, the panels of Fig. 8 together convey this prior dependent knowledge of DA from our quasi-DA matrix element.
We repeated the above fits for all analysis choices, which now includes the truncation order N G = 2, 3 and 4 in Eq. (46). In the top panel of Fig. 9, we show our estimate of the pion DA as a function of u, after taking into account all the analysis variations, and summarize them with the statistical and systematic error bands. To be cautious, we present the reconstruction using a relatively broad prior width δ = 0.2 on the expansion coefficients. Nevertheless, the reconstruction in the case of DA is sensitive to the value of δ, however large it is, and hence, one should interpret the reconstruction of DA in Fig. 9 as a specific u-dependence, given a somewhat broad prior. We compare our result with the asymptotic DA shown as the black dashed curve. Within the precision allowed at δ = 0.2, we can only resolve an overall flat DA over a range of u ∈ [0.2, 0.8] with sharp fall offs, u α and (1−u) α with α ≈ 0.3, to 0 on either side. If one focuses only on the central value of the reconstructed DA, one sees a tendency for a platykurtic DA as noted in Refs [97][98][99]. The lattice data does not have the sensitivity to further resolve the concavity or convexity within the flatter regions, unless one is willing to impose a more stringent prior width δ. Apart from providing a reconstruction of the DA, the ansatz based analysis also provides a way to estimate the moments of DA. The usage of ansatz can be thought of as a way to regulate the values of moments at larger-n for which the lattice data is less constraining, and therefore, provides robust values for smaller-n moments. From the Ansatz based analysis above with δ = 0.2, we estimate the Mellin moments as By comparing the values with Eq. (42), we see that the ansatz based reconstruction for x 2 agrees quite well with the completely model-independent reconstruction. The estimates of x 4 also agree with each other, however, the usage of ansatz has substantially reduced the error. Thus, from both the model-independent moments analysis and the model-dependent reconstruction analysis, we find the values of x 2 and x 4 to be the quantities that we could reliably extract from our lattice data.
As another way to summarize our results with less modeling artifacts, we present the MS light-front ITD corresponding to the pion DA in the bottom panel of Fig. 9 in the range of λ that we have lattice data for and performed our analysis on. To infer the MS ITD, we used Eq. (11). Since we need only the information on the Mellin moments to construct the light-front ITD, we show the resultant ITD based on the above Ansatz-based analysis as the red band enclosed between the dot-dashed lines, and the result based on the Mellin moments analysis in the previous subsection as the blue band enclosed between solid lines. In both cases, the darker inner bands are the statistical error bands whereas the lighter outer bands include both statistical and systematical errors. We see that both the model-independent and ansatzdependent reconstructions have similar behavior in the range of λ that is constrained by the lattice data, with the latter being a more precise determination. We show the light-front ITD corresponding the asymptotic DA as the black dot-dashed curve. We also compare our result with the expected ITD for a flat pion DA shown as the magenta dot-dashed curve. Our result is clearly below the asymptotic DA expectation, and closer to the expectation from a flatter DA. This expectation is a less model-dependent manifestation of the DA reconstruction seen in the top panel. The key quantities that characterize exclusive QCD processes, such as such as the photon-pion transition form factor [10][11][12][13], electromagnetic form factors and GPDs can be factorized into convolutions of DAs and the perturbatively-calculable partonic hard-scattering amplitudes if the momentum transfer is sufficiently large. For electromagnetic form factor this factorization was introduced long time ago [5,7,25]. For photon-pion transition it was discussed in Refs. [14,15], and for gravitational form factors it was discussed in Refs. [100,101], while for GPDs it was discussed in Refs. [102,103]. The energy scale where this leading-twist DA-based factorization may work is unknown at present. This an important question that can only be answered by experiments or through lattice QCD computations. In this subsection, we put aside this question and simply make predictions for electromagnetic and gravitational form factors of the pion based on leading twist factorization and our DA results. These predictions can be compared to the lattice or experimental results at large momentum transfer and clarify the range of applicability of the leading twist factorization for the form-factors.
At large Q 2 , the pion electromagnetic form factors F π (Q 2 ) can be factorized as where Q 2 is the momentum transfer and T F is the hardprocess kernel. Though the form factor F π (Q 2 ) is scale independent, the fixed-order perturbative factorization introduces dependence on both renormalization and fac-torization scales µ 2 R and µ 2 F . Here Φ(u, µ 2 F ) is defined as, where f π is the pion decay constant discussed in Section V A. At leading order (LO), the hard kernel reads [106], withū = 1 − u and the running coupling constant where β 0 = 11 − 2 3 n f , and we use n f = 3, Λ QCD = 0.2 GeV in this paper. We take our model fit result at the initial scale µ 0 = 2 GeV, and evolved it to µ F by first expanding φ(u, µ 0 ) in the Gegenbauer basis in Eq. (45) up to a sufficiently large order n = n max = 20, and then evolving those Gegenbauer moments from a n (µ 0 ) to a n (µ F ) using We choose µ R = µ F = Q as the central value of the scale setup, and vary the renormalization scale µ R by a factor of 2 to estimate the perturbation uncertainty. For the ease of implementation, we used the 1parameter Ansätze φ 1−param (x, µ R ) from our analysis using (N HT , N LC , n 0 3 , z min 3 , z max 3 ) = (1, 1, 1, 2a, 8a). The results are shown in Fig. 10 with statistical error bands and compared with the experimental data from the F π collaboration [22] as well as the calculations from the Dyson-Schwinger equation (DSE) [104] and Minkowskispace Bethe-Salpeter equation (BSE) [105]. As one can see, our prediction using the LO kernel is systematically lower than the DSE and BSE calculations. It was found that the matched form factors could increase with NLO corrections [106], and the higher-twist corrections may also make a significant contribution [107]. However, all these arguments can only be tested by the future experimental results with large momentum transfer Q 2 up to 6 GeV 2 of the JLAB E12-09-001 experiment [24] and up to 40 GeV 2 of the new Electron-Ion Collider (EIC) facility [2].
The gravitational form factors (GFFs) of the pion are the transition matrix elements of the QCD energy momentum tensor, Though T µν QCD is conserved, and is therefore UV finite and scale independent, the quark and gluon contributions, T µν q and T µν g , are not and depend on the renormalization scale µ R . This dependence is governed by the   11. The gluon GFFs Q 2 Ag(t), quark GFFs Q 2 Aq(t) and Q 2 Cq(t) as well as the pion trace anomaly P | β(g) 2g F 2 |P at large −t predicted by our determination of pion DA are shown. The bands come from the statistic errors and we vary the renormalization scale µR by a factor of 2 to estimate the perturbation uncertainty. The direct lattice calculation of pion gluon GFF Ag(−t) [108] (the multipole fit result) using unphysical pion mass mπ = 450 MeV is shown (MIT21, the yellow band in the top-left panel) for comparison.
corresponding anomalous dimension. The gluon gravitational form factors for the pion can be parametrized as, whereP = (P + P )/2 is the average momentum, ∆ = P − P is the momentum transfer, and −t = −∆ 2 = Q 2 . Similar to the case of the electromagnetic form factor, the leading-twist GFFs perturbative factorization reads, where µ F -dependence can be introduced in A π g due to the use of a fixed-order hard kernel. At leading order, the kernels are [100], Due to the traceless feature of Eq. (55) one also has the relation C π g = −C π u+d = t 4m 2 C π g . In terms of the same factorization formula, the hard coefficients A q and C q in the quark sector are, In addition, the gluon scalar FF defined as, is closely related to the trace anomaly [109], Using the same factorization convention as Eq. (56), the leading-order hard kernel reads [101], Comparing to the tensor GFFs A q,g and C q,g shown above, one can observe that G π q does not have a 1/(−t) pre-factor and therefore will become flat for large −t. In Fig. 11, we show the perturbatively determined GFFs and pion trace anomaly at large −t as expected from our determination of pion DA. The bands come from the statistical errors and we vary the renormalization scale µ R by a factor of 2 to estimate the perturbation uncertainty. The direct lattice calculation of the pion gluon GFF A g (−t) [108] (the multipole fit result) using unphysical pion mass m π = 450 MeV is shown (MIT21, the yellow band in the top-left panel) for comparison. And it can be seen our estimate of the perturbative contribution is much smaller, which could come from sizable higherorder perturbative or higher-twist contribution. Direct lattice calculations or experimental results at large −t are needed to clarify the issue.
In order to perform the above perturbative convolutions, we relied on the Ansatz-based reconstruction to a large extent. Instead, one could ask if there is a way to perform an alternate less model dependent analysis based only on the ITD in the range of λ spanned by the lattice data, or equivalently, only using the moments that the lattice data is sensitive to. To address this question, we can refer to the LO perturbative convolution in Eq. (49) for electromagnetic form-factor that makes use of the integral, Using Eq. (45), one can see that, u −1 = 3 ∞ n=0 a 2n (µ). If one truncates the sum over Gegenbauer moments up to n = 1 at µ = 2 GeV, then one finds u −1 = 3.72 (45), whereas when one sums over 20 Gegenbauer moments using the 1-parameter Ansätze, we get u −1 = 4.87 (19). As an alternative, we can write Eq. (63) using only the ITD as, With λ max ≈ 6 as in this work (see bottom panel of Fig.  9), we find this region of I(λ, µ) to contribute 2.64(2) to u −1 , which is only about 50% to the Ansatz-based expectation for u −1 , and the rest comes from λ > λ max . Therefore, for the two model-independent methods to be reliable, one has to truncate at much higher Gegenbauer moments or larger λ, which poses a challenge for lattice calculations. As a result, one has to rely on the Ansätze for the x-dependence of DA, but the systematic uncertainty from truncation is transformed to the model dependence of the Ansätze. It would be important in the future to compare and cross-check our current results for the form factors here with the expectations based on xspace LaMET DA matching on the same ensemble.

VII. CONCLUSIONS
We presented a lattice QCD study of the quasi-DA matrix element in real-space using the leading-twist OPE method for the first time. We performed our study at the physical point using clover-improved Wilson valence quark propagators determined on a physical HISQ ensemble. The quantities central to this work are the renormalization group invariant ratios of quasi-DA matrix elements, with non-zero momenta in both the numerator and denominator; the non-zero momenta were inevitable as well as helped us remain closer to the leading-twist approximation. In the first part of the paper, we adapted the results from Refs. [43,75] and presented the analytical perturbative results for the leading-twist expansion in the forms of the Mellin OPE at next-to-leading order and the conformal OPE at leading-log order. These expressions formed the basis for our determination of the pion DA from quasi-DA matrix elements.
From the leading-twist description of z · P and z 2 dependencies of the ratios of quasi-DA matrix elements, we extracted the Mellin moments and captured the xdependence of the pion DA based on fits to various Ansätze. We first checked the validity of leading-twist dominance in our matrix element using a fixed-z 2 analysis. Then, from a model-independent determination via fits to the few lowest Mellin moments at a factorization scale µ = 2 GeV, we were able to obtain a relatively precise determination of the second Mellin moment, x 2 = 0.287 (6)(6), and for the fourth Mellin moment, we obtained x 4 = 0.14(3)(3); the first parenthesis gives the statistical error and the second one specifies the systematical error coming from variations in various analysis choices, such as fit ranges for z 2 . Based on NLO perturbative evolution of our corresponding results at different µ evolved to µ = 2 GeV, we estimated our perturbative uncertainty to be within our combined statistical and systematical errors. We reached a similar conclusion from the differences in the estimates of the Mellin moments from the Mellin-OPE and Conformal-OPE. We found the Ansatz based reconstruction of the x-dependence of the pion DA to be sensitive to the model used; by using a complete set of functions to expand the DA and by imposing constraints on their expansion coefficients, we systematically reconstructed the pion DA. Using a weak constraint, we found the DA at µ = 2 GeV to be flatter than the asymptotic DA in the region around x = 0 (or equivalently, u = 0.5). From the Ansätze-based reconstruction of the pion DA, we found the expected the large Q 2 dependence of electromagnetic and gravitational form factors using the leading-twist LO convolutions. It would be interesting in the future to compare the values of the pion form-factors at these large Q 2 with the perturbative expectations.
The systematical error in this work stemmed only from the choices of fit ranges, type of higher-twist corrections and other such analysis choices. Another source of systematic error could be due to finite lattice spacing corrections. Since, we used an ensemble at a fixed lattice spacing, a = 0.076 fm, we were unable to quantify the effect in this work, and we need to revisit this in a future work. We found the perturbative uncertainties to be about 3% as estimated through differences in results for x from Mellin-and Conformal-OPE, and through the effect of evolution to 2 GeV starting from different initial scales used for fits. However, in the present work, we were not able to directly address this issue using NNLO DA matching as is the state-of-art for the current lattice PDF calculations. In the immediate future, we plan to extend the current work using the leading-twist expansion of the quasi-DA to study the Kaon DA and quantify the effects of explicit SU(3) flavor symmetry breaking. In this appendix, we discuss the calculation of the renormalization of the vector current Z V and axial-vector current Z A in RI-MOM scheme for our setup with a = 0.076 fm. We use off-shell quark states in the Landau gauge with different values of lattice momenta ap µ = 2π L µ (n µ + 1 2 δ µ,0 ).
To minimize the discretization errors the lattice momenta are substituted by ap µ = sin(ap µ ), so the renormalization point is (ap R ) 2 = µ=1,4 (ap µ ) 2 . In Fig. 12 results for Z V (p R ). The vector current renormalization constant should not depend on p R , because in the a → 0 limit the local current is conserved. Nevertheless, we see a significant dependence on p R . This dependence is caused by non-perturbative effects, that for large values of p R can be parameterized by local condensates. As we use off-shell quark states in the Landau gauge in the RI-MOM renormalization procedure, the lowest dimension local condensate is the dimension-two gluon condensate A 2 [110,111]. Lattice artifact shows up as the breaking of the rotational symmetry on the lattice. We see from Fig. 12 that the fish-bone structure in the lattice data at the level much larger than the statistical errors on Z V for large values of ap R . Therefore, to obtain Z V we fit our lattice data with the following form: where ∆ (4) = µ (p µ ) 4 p 4
(A3) This form is motivated by the 1-loop lattice perturbation theory [112,113] and the perturbative analysis with dimension two gluon condensate [114]. For the nonperturbative clover action k = 2, while for Wilson action k = 1. For HYP smeared clover action with tadpole improved value of c sw we expect O(a) discretization errors to be proportional to α 2 s with a very small coefficient, so it is reasonable to assume that the dominant cutoff effects scale like a 2 . Nevertheless we also perform fits using k = 1. To limit the size of the lattice artifacts we impose the additional constraint: ∆ (4) < 0.4. We performed different fits varying the fit interval in p R as well as setting some coefficients to zero in certain cases. Fits with C = 0 typically have very large χ 2 but this has almost no effect of the extracted Z V value. From the fits we obtain Z V = 0.947 (8), where the error is mostly systematic and corresponds to the scattering of the results from different fits. We could also estimate Z V from the matrix element of the vector charge of the pion calculated in Ref. [85]. Using the result for the matrix element from the two state fit we obtain Z V = 0.9534(5) [85]. This agrees with the above result within errors.
In Fig. 12 bottom panel, we also show the ratio Z A /Z V as function of p R . This ratio too should be independent of p R . We see some dependence on p R due to non-perturbative effects, though it is considerably milder than for Z V . This is likely due to the fact that the leading non-perturbative contributions cancel out in the ratio Z A /Z V . The lattice discretization effects also seem to largely cancel in the ratio Z A /Z V and no clear fishbone structure can be seen in our data. The p R dependence of Z A /Z V is incompatible with B/p 2 R form. Therefore, we fit our data with const + B /(ap R ) 4 form. This gives Z A /Z V = 1.0168. For (ap R ) 2 > 1.5 it is also possible to fit the data with constant, which gives Z A /Z V = 1.01514. Combining these results with the value of Z V from the pion matrix element of the vector charge we obtain Z A = 0.969 (1). The error is systematic and is due to the difference of the two fits of Z A /Z V .

Appendix B: Details on the Mellin Moments fit
In Section VI C, we presented the results on the fitted values of the first few Mellin moments based on fits to leading twist Mellin OPE along with lattice correction of the form in Eq. (39) with l 2 being a fit parameter, and with higher twist corrections in Eq. (40) with h 0 and h 1 as possible fit parameters. In this appendix we discuss the results for these additional fit parameters in the Mellin OPE fits with no positivity constraints on the moments. In the first three panels of Fig. 13, we show the scatter of best fit values for l 2 , h 0 and h 1 respectively for various fit types specified by (N HT , N LC , n 0 3 , z min 3 /a, z max 3 /a) beside the data points. Since certain parameters are not part of all fit types (e.g., h 1 does not occur when N HT = 0, 1), data points for those cases are left missing in the different panels. From the figure, the conclusions specified in the main text become clearer; namely, the presence of l 2 has only a marginal effect, whereas the effect of the highertwist term h 0 cannot be neglected. Also, N HT = 1 is sufficient for the fits for the present data quality, whereas inclusion of h 1 with N HT = 2 only makes the fits noisier. In the rightmost panel of Fig. 13, we show the scatter of minimum χ 2 /df in the various fits. The error bars on the minimum χ 2 /df is obtained from the Jackknife blocks. The mean values of χ 2 /df for different fits all lie approximately between 1 and 1.6.