Precision Control in Lattice Calculation of $x$-dependent Pion Distribution Amplitude

We present a new Bjorken $x$-dependence analysis of a previous lattice quantum chromodynamics data for the pion distribution amplitude from MILC configurations with three lattice spacing $a=0.06,0.09, 0.12$~fm. A leading renormalon resummation in renormalization as well as the perturbative matching kernel in the framework of large momentum expansion generates the power accuracy of the matching to the light-cone amplitude. Meanwhile, a small momentum log resummation is implemented for both the quark momentum $xP_z$ and the antiquark momentum $(1-x)P_z$ inside a meson of boost momentum $P_z$ up to 1.72 GeV along the $z$ direction, allowing us to have more accurate determination of the $x$-dependence in the middle range. Finally, we use the complementarity between the short-distance factorization and the large momentum expansion to constrain the endpoint regions $x\sim 0, 1$, thus obtaining the full-range $x$-dependence of the amplitude.


Introduction
Distribution amplitudes (DAs) are important observables for both theoretical and phenomenological reasons within the realm of quantum chromodynamics (QCD). The DA of a meson describes the probability amplitude of identifying the meson in a quark-antiquark Fock state on the lightcone, carrying longitudinal momentum fractions x and 1 − x, respectively. It is also known as the leading Fock wave function of the meson. They are important as inputs to many exclusive processes with large momentum transfer, such as the B-meson decay, that can be factorized into the nonperturbative DA and the hard-scattering kernel [1,2]. Although the DAs are important quantities in QCD, their properties, such as the moments, the shape and the endpoint power-law behavior are still undetermined from experiments [3,4,5,6]. A direct nonperturbative calculation of the DAs from lattice QCD is thus of great interest.
The nonperturbative physics of partons is defined on the lightcone, i.e., in the effective limit of infinite momentum. Direct calculations on the lightcone are inaccessible on the lattice due to the dependence on real time. Early calculations determined DAs by calculating their lowest moments from local twist-2 operators [7,8,9,10,11,12] or from nonlocal current-current and quark bilinear correlators [13,14,15,16,17]. The local-operator calculations provide precise measurements up to the second moment of the DA [12], but the increasing noise and the nontrivial mixing in the lattice renormalization make it very difficult to access higher moments. The nonlocal-operator calculations analyze data in a certain current-current displacement or Wilson-line length z range, where the short-distance factorization is valid; this either allows us to obtain the lowest few moments, or needs a model assumption to fit the x-dependence [15,17]. A direct x-dependence calculation has not been possible in the two traditional methods.
The method of large-momentum effective theory (LaMET) [18,19,20] offers a different approach, which starts from the Euclidean matrix element of equal-time, spatially separated fermion fields. After renormalization, we can physically extrapolate these matrix elements to large distances and Fourier transform them to momentum space. We use field-theoretical large momentum expansion to match the data at finite hadron-momentum to the light-cone lightcone DA, but not the local x-dependence. The two methods complement each other [34], enabling us to combine the local information from the LaMET calculation and the global information from the short-distance OPE. We model the x-dependence outside the region x ∈ [x min , x max ] and fit to the short-distance correlations, obtaining a modelindependent mid-x distribution and a model-dependent endpoint distribution. The endpoint distribution is constrained by the physical requirement ϕ(x) → 0 when x → 0 or 1 as well as the requirement of continuity with the LaMET calculation, which limit the model dependence. This combined approach provides a full x-dependence calculation.
With the above three ideas (elimination of the linear correction, resummation of large logarithms and constraining the endpoint regions), we improve the analysis of the lattice quasi-DA data to extract the full x-dependence of the lightcone DA with improved accuracy. The rest of the article is arranged as follows. In Sec. 2, we describe the DA calculation using LaMET, discuss the ambiguities in the renormalization and perturbative matching, and present how we achieve power accuracy in the LaMET matching. In Sec. 3, we discuss the origin of two different physical scales in the LaMET matching and show how to resum them. In Sec. 4, we discuss how to use the short-distance OPE to constrain the endpoint regions to extend our calculation to the full range of x. In Sec. 5, we apply LRR renormalization and LRR matching with RGR to extract the lightcone DA, then use complementarity to obtain the full x dependence. Finally, we conclude in Sec. 6.

Renormalization and Power Accuracy
The correlator that defines the pion DA on the lightcone is where W(0, η − ) =P exp −ig η − 0 ds n µ A µ (ns) is the Wilson line between the two points 0 and η − , andP is the pathordering operator. Lightcone coordinates are defined for a general Lorentz vector V µ , as V ± = 1 √ 2 (V 0 ± V 3 ), since we may assume without loss of generality that the meson is traveling in the z(≡ η 3 ) direction. The term g is the coupling, A µ denotes the gauge field and |π(P)⟩ is a pion state with 4-momentum P + . The variable x is the fraction of the meson momentum carried by the constituent parton.
The operator in Eq. (1) has dependence on real time and is, thus, inaccessible directly on the lattice. The method of LaMET begins with the following "quasi" correlation: whereh R π is the renormalized coordinate-space matrix element defined in the second line. The lightcone DA, ϕ π (x, µ), is related to the quasi-DA (qDA) in the large momentum P z limit viã where the C(x, y, µ, P z ) is the perturbative matching kernel in momentum space, and the residual quadratic in Λ QCD /P z comes from higher-twist effects and is only leading if the linear divergence in the bare operator in Eq. (2) were not present as we explain below. The bare matrix elements we compute on the lattice are theh B π (z, P z , a) terms corresponding to Eq. (2) before renormalization, so our data are initially in position space and contain UV divergences. The spatial Wilson line W(0, z) has a linearly-divergent self energy of size 1 a , so besides the usual logarithmic divergence, the linear divergence must also be removed through a multiplicative renormalization [35] before extrapolating to the continuum.
where Z R (z, a) ∼ e −δm(a)z is the renormalization constant with the linearly divergent mass counterterm δm(a) ∼ 1 a . When renormalizing the linear divergence, one could in principle also choose to subtract a finite constant term along with it. The choice of this finite piece defines the renormalization scheme. Also, when expanding δm(a) as a perturbation series in the strong coupling, α s , δm = 1 a n r n α n+1 s (a −1 ), the coefficient r n ∼ n! grows factorially at higher orders due to an infrared renormalon effect [36,37]. Thus, the series is divergent for any α s , and the sum is ill-defined. To fix this degree of freedom, we need to introduce an additional renormalization scheme for the linear divergence to define δm(a, τ) unambiguously, with a new τ-dependence. This result varies by O(Λ QCD ) in different τ-schemes, so an ambiguity of O(zΛ QCD ) arises in the renormalization factor. The same intrinsic ambiguity appears when we try to extract δm from fitting lattice data, where δm is always mixed with another non-perturbative quantity, such that we have a freedom to choose among different fitting results. Similarly, a calculation of the perturbative matching kernel C(x, y, µ, P z ) also suggests a factorial growth with the same pattern [38]. The lightcone distribution ϕ(x, µ) is obtained by convoluting the inverse matching kernel C −1 (x, y, µ, P z ) with the renormalized quasi-DA,φ π (x, P z ), both containing the ambiguities. The combination will, in general, result in a linear correction O Λ QCD xP z to the matching [29], where we have ignored the (1 − x)P z scale for simplicity, which can be recovered by a substitution x ↔ 1 − x, due to the symmetry of the matching. When the hadron momentum is large enough, this correction is not important. But the hadron states in lattice calculations are usually moving with P z ∼ GeV, where the linear correction can be large, especially near the endpoints, and more important than the quadratic higher-twist effects. In principle, these ambiguities from the renormalization and the perturbative matching can cancel because the twist-2 lightcone DA ϕ(x, µ) is free of the linear divergence and the infrared renormalon. Thus, the power accuracy up to O Λ QCD xP z is the best we can achieve without knowing higher-twist information. It is shown that they indeed cancel only when the renormalization of linear divergence and the regularization of the matching coefficients are defined in the same τ-scheme [32]. To achieve this accuracy, we need to carefully define the renormalization scheme for the linear divergence δm and regularize the perturbative matching consistently to eliminate the linear correction.
A recent work [25] uses a fixed-order approach to handle this ambiguity by introducing an additional twist-three mass parameter, denoted by m 0 , in the renormalization process to ensure that the short-distance behavior of the renormalized matrix element is in agreement with perturbation theory. The approach is still not good enough for several reasons: 1) A bridge is missing to connect the lattice calculation and the perturbative calculation, usually known as the scheme-conversion factor in the renormalization; 2) Resumming the logarithms ln z 2 µ 2 at short distances clearly suggests that m 0 is not a constant but has a large dependence on z, mainly due to the fixed-order truncation not being a proper scheme to regularize the divergent series [32]; 3) A fixed-order matching is used, which cannot eliminate the linear ambiguity.
We propose a new approach, aimed at eliminating such a correction to achieve the power accuracy, as demonstrated in Ref. [32]. It includes four steps: • Modify the perturbative matching coefficients through a leading renormalon resummation (LRR) with a principal value (PV) prescription, defined as the τ-scheme.
• Determining the non-perturbative twist-3 parameter m 0 (τ) through the matching condition that the renormalized P z = 0 lattice data agree with the LRR-improved Wilson coefficients in the τ-scheme up to twist-3 accuracy in the OPE at short distances; • Renormalize the P z > 0 lattice data with the m 0 (τ) extracted from the previous step in the τ-scheme; • Extract the DA with the LRR-improved perturbative matching kernel.
The m 0 (τ) parameter is fixed by P z = 0 pion quasi-PDF data and used for the renormalization of the quasi-DA at nonzero momentum. The justification for this choice is that the linear correction from the ambiguity in the linear divergence is independent of the momentum of the external state and the Dirac structure. δm is universal for the Wilson-line self energy, which is the same in these observables (P z = 0 and P z > 0) obtained from the same gauge action; the leading renormalon contribution that we resum in the MS perturbative calculation also originates from the Wilson line self energy, thus is the same for these observables, up to an overall phase factor depending on the external states' momentum. Thus, once the cancellation of the ambiguities is achieved for one observable, it is also guaranteed for other observables of a similar structure, i.e., with the same Wilson line in the quark bilinear operator but with a different Dirac structure or different external states.
The LRR improves the renormalization method in two aspects, as we will show in Sec. 5. Firstly, the e m 0 (τ)z factor is extracted with LRR, so it is different from that extracted in fixed-order perturbation theory. With LRR improvement, its extraction is almost independent of small-z values, and is determined with a significantly reduced uncertainty from scale variation [32]. Secondly, in the hybrid scheme, the renormalized matrix elements are divided by P z = 0 perturbative results (i.e., the Wilson coefficient C 00 ) at short distances, whose z-dependence is improved after LRR in the sense that they are more consistent with the OPE at short distances, with the lowest few moments as inputs. This moment is supposed to be consistent with the one extracted from a renormalization-independent ratio between two different momenta, which will be discussed in more detail in Sec. 4. Thus, a comparison between the moments extracted from the renormalization-dependent matrix element and from the renormalization-independent ratio will test whether the renormalization is properly done.
The idea of LRR is to resum the leading factorially divergent high-order terms to all orders in the perturbation series. Then the remaining part of the series, if without other renormalons, is convergent, and any leading power correction could be fixed by a regularization of the resummed divergent part. Although it is impossible to analytically calculate the perturbation for specific processes to all orders, we can calculate a specific type of bubblechain diagrams [39] in the large β 0 limit. Beyond the large-β 0 limit, the asymptotic form of the leading renormalon pole is known [36,40], whose overall strength has been estimated from perturbation series of the heavy quark pole mass [40,41] and lattice calculations of the static potential [37]. Thus we can also choose to resum these known asymptotic forms. These two approaches both resum the leading pole corresponding to the linear divergence, but have different "background" effects that are higher powers of Λ QCD and higher order of α s . Thus they are supposed to make slightly different predictions in the mid-x region, as we will discuss in Sec. 5.

LRR in the large β 0 limit
For quasi-PDF operators, a calculation for bubble-chain diagrams has been done in Ref. [38]. Note that only the Wilson-line self-energy diagram (also called the "tadpole" diagram in Ref. [42]) is relevant to the leading renormalon, so we can ignore the other diagrams which only account for higher renormalon poles.
By resumming the tadpole diagrams, the LRR in the large-β 0 limit modifies the P z = 0 matrix element for the DA, i.e., the Wilson coefficient C 00 (z, µ), in the following way: where C tp 00 (z, µ)| PV is the resummed diagrams with the principal value prescription for the poles defined as scheme τ, and C tp,(i) and the corresponding Wilson coefficient [42] C (1) 00 (z, µ) = α s C F 2π The P z > 0 matrix elementsH(z, P z , µ) are corrected by LRR in a similar way, where the momentum dependence only enters through a phase factor, which can be Fourier transformed to obtain the correction to the NLO matching kernel in MS, where the first part is not a traditional convergent function, but a distribution operating on the DA function through a convolution, whose effect is convergent mathematically. In practice, it is enough to perform a truncated numerical evaluation to some large z max , e.g., z max = 10 fm.
In the ratio scheme [43,44,45], the ratio between two momentums is free of linear divergence, thus no LRR modification is needed, and the matching kernel C ratio is unchanged. In the hybrid scheme, the correction is an integration from z s to z max during the Fourier transformation, where the first term can be calculated numerically to z max in practice. It is also straightforward to derive the LRR correction to the DA Wilson coefficients by expanding Eq. (11) in zP z : which can be applied to the OPE of short distance correlations.

LRR of the asymptotic series
Besides resumming the leading renormalon pole, the large-β 0 approximation introduces extra effects in subleading renormalon poles. Alternatively, as discussed in Ref. [32], we can resum the asymptotic form of the leading renormalon contribution, which only includes the leading renormalon pole. In this approach, we utilize the fact that the leading renormalon contribution originates from the heavy quark pole mass m = µ n r n α n+1 s , with a known asymptotic form in large perturbation order n [36,41,37], where b 0 = β 1 /2β 2 0 and c 1 = (β 2 1 − β 0 β 2 )/(4b 0 β 4 0 ) are from higher orders in the QCD beta function. Using an analytical method in Ref. [40], the overall strength can be determined as N m (n f = 3) = 0.575, N m (n f = 4) = 0.552. Thus the contribution to the DA Wilson coefficients has the following form at large n: Similar to the LRR in the large-β 0 limit, we can resum the asymptotic form with the PV prescription, It's easy to verify that the ambiguity of this integral is linear in zΛ QCD and independent of µ. Note that the Fourier transformation of this correction can be calculated analytically, but the explicit linear-z dependence will be transformed into a singular distribution of x − y, including derivatives of the δ(x − y) function. It is numerically very unstable if this function is applied to discrete data. So a regularization is applied, by multiplying the linear-z term with a small exponential decaying factor exp(−ϵ m z), which will result in extra higher-twist corrections O With such a regularization, the correction to the hybrid-scheme matching will be where ∆ xy = |x−y|. The overall factor C asympt kl (z, µ) PV /z−r 0 µα s (µ) only depends on µ and can be integrated numerically. The total correction is written as a plus function to guarantee the current conservation because one term proportional to δ(x − y) has been omitted. Testing with some different ϵ m ∼ 20 − 100 MeV values, and with ∆x = 0.01 as the step size of our numerical methods in the momentum space matching, we find the results are consistent and stable. Working with smaller ϵ m requires a finer discretization of the data as a function of x or y.

Resummation in coordinate space
To study the resummation of large logarithms, we start from a simpler case, the coordinate-space matching of the quasi-DA. It is more straightforward because only one physical scale, z = λ/P z , is involved in the matching. The renormalized quasi-DA matrix element,h R (λ, P z ), can be matched to lightcone DA, h(λ, µ), through where Z(ν, z 2 , µ 2 , λ) is the perturbative matching kernel in coordinate space. In the ratio scheme [43,44,45], where L = ln z 2 µ 2 e 2γ E /4 is the only scale-dependent logarithm appearing in the kernel and C F is the quadratic Casimir for the fundamental representation of SU(3). At either short distances, z → 0, or long distances, z ≫ Λ −1 QCD , the logarithm becomes large. We can eliminate the logarithm by setting µ 0 = 2e −γ E z −1 on the right-hand side of the equation. Then an RG evolution to the default scale, e.g. µ = 2 GeV, will resum the large logarithms at that scale, whereV is the coordinate space representation of the Efremov-Radyushkin-Brodsky-Lepage (ERBL) evolution kernel [46,47,48,49], Both the evolution and matching kernels can be made purely real by multiplying by a phase factor e iλ(1−ν)/2 , and applying these to the phase-rotated matrix elements which is purely real for symmetric DAs that satisfy ϕ(x) = ϕ(1 − x). So the matching and the evolution preserve the symmetry of the DAs. Such a resummation works fine at short distances, but at long distances the scale µ 0 hits the Landau pole, indicating that the perturbation theory as well as the entire short-distance operator expansion break down. Without knowing the correct information at large z, we are unable to extract the x dependence of the DA.

Origin of two different scales
In large momentum expansion, we perform the resummation and matching in momentum space, and the highertwist non-perturbative physics appear now at the endpoint regions x → 0 and x → 1, which we will choose to model using complementarity, as we will discuss in the next section. In momentum space, the quasi-DA is matched to the lightcone DA throughφ The momentum-space matching kernel C(x, y, µ, P z ) can be obtained from a double Fourier transformation of the coordinate-space matching in Eq. (19), To trace how the physical scale and the logarithm transform, we can check the double Fourier transform of the logs, L, in Z(ν, λ 2 y 2 P 2 z , µ 2 ) of Eq. (19). The terms involved include f (ν)L and f (ν)Le −iλ(1−ν) . In dimensional regularization d = 4 − 2ϵ, higher-order logs L n can be expressed as the O(1) term of n! ϵ n µ 2 z 2 e 2γ E /4 ϵ in the ϵ expansion. Integrating λ first, we obtain which can be expanded in ϵ → 0 [42]. Note that only when 0 < x/y < 1 is it possible for |x/y − ν| to be zero in the integration region ν ∈ [0, 1], and then the expansion of |x/y − ν| −1−2ϵ generates the leading divergent term −1 2ϵ δ(x/y − ν). So when x is in the nonphysical region, or when 1 > x > y > 0, the expansion does not contain any leading logarithms of (ln µ) n . When 0 < x < y < 1, the expansion yields additional log terms 1 0 dν f (ν) ln n (µ 2 ) − n ln n−1 (µ 2 ) ln 4y 2 The remaining integral preserves the structure of the log and only changes its coefficients. Thus, we get the physical scale 2yP z for this term. On the other hand, the other term f (ν)Le −iλ(1−ν) after a double Fourier transformation becomes with x ≡ 1 − x and y ≡ 1 − y, which does not contain any leading logarithm when 1 − x is nonphysical, or when 0 < x < y < 1. When 1 > x > y > 0, the ϵ expansion yields So the two different physical scales correspond to different regions of x and y. The scale 2yP z for x < y corresponds to the quark-splitting process; the scale 2(1 − y)P z for x > y corresponds to the antiquark-splitting process. Note that the two scales 2yP z and 2(1 − y)P z we obtained at the current stage both depend on the convolution variable y, which is, in principle, not implementable because we cannot have different scales µ for different y in ϕ(µ, y). However, note that the matching kernel C(x, y, µ, P z ) is almost localized, i.e., the region of x ∼ y is greatly enhanced compared to any other regions. As a result, the proper scale choice to resum the RG logarithm would be 2xP z and 2(1−x)P z instead, corresponding to the quark and antiquark momentum fractions in the quasi-DA. Moreover, a comparison between the quasi-PDF and DIS shown in Ref. [33] suggests that 2xP z is the proper scale in the quasi-PDF case, also supports the scales to be 2xP z and 2(1 − x)P z in our quasi-DA. We can examine the sensitivity to this choice of scale by slightly varying the value from 2xP z to 2cxP z with c ∈ [0.75, 1.5], which roughly correspond to a ±30% change in α s near 1 GeV.

Resummation of two different logarithms
The perturbative matching kernel in Eq. (23) C(x, y, µ, P z ) = δ(x − y) + C (1) (x, y, µ, P z ) + O(α 2 s ) has been calculated to 1-loop order [50], where C (1) is the O(α s ) term of the matching, and satisfies the quark-antiquark symmetry Moreover, it contains logarithms of both kinds as discussed in the previous subsection, which becomes large at the end-point regions. Indeed, in the two regions x < y and y < x, the matching kernel has different µ dependencies, corresponding to the piecewise function of the ERBL evolution kernel: where and the plus function is The logarithms become large in the matching kernel C(x, y, µ, P z ) for x close to both endpoints x → 0 or 1, so a resummation of large logs is necessary. The traditional method of resummation is to choose a scale µ 0 in the scale-independent factorization such that the large logs in C(x, y, µ 0 , P z ) are eliminated, and apply the RG evolution of ϕ(y, µ): where P is the path ordering of the evolution path, andV is the operator corresponding to the ERBL kernel acting on ϕ. However, such an approach does not work in the quasi-DA's case, because there are two different scales in the piecewise matching kernel. The logarithm in x < y is ln[4x 2 P 2 z /µ 2 ], corresponding to the quark-splitting; the logarithm in x > y is ln[4(1 − x) 2 P 2 z /µ 2 ], corresponding to the antiquark-splitting. No single choice of µ 0 is able to eliminate the two logs at the same time. We have to develop a different strategy to resum the log in the DA.
We start from separating the matching formula in two different regions ϕ(x) = ϕ(x, µ) + x<y dy C (1) (x, y, µ)ϕ(y, µ) + x>y dy C (1) where we label the integral in two regions as C (1) L and C (1) R convoluted with ϕ. One idea is to set the two terms to different scales, C (1) L (2xP z ) ⊗ ϕ(2xP z ) and C (1) To make this possible, we need to split the first term into two parts and combine with them separately. After the split, the two parts need to be scale-invariant individually up to order α s . So we need to find the weights, w L (x) and w R (x) = 1 − w L (x), for the following split such that the two lines are individually scale-independent. If we force w L (x) to be scale independent, then simple algebra gives where we define V (0) L ≡ V (0) | x<y and V (0) R ≡ V (0) | x>y in a similar way. The solution w L (x) is scale-dependent at order α s , contradicting our requirement. However, noticing that ϕ(µ) andφ differ by O(α s ) corrections, and the quasi-DAφ is scale independent, we may substitute the ϕ(µ) in Eq. (36) withφ, making w L (x) scale independent. The new solution satisfies the scale dependence of the original equation at order O(α s ), which cancels the µ dependence in C (1) L (µ) ⊗ ϕ(µ). So we write the matching formula as where the two scales are chosen to be µ 1 = 2xP z and µ 2 = 2(1 − x)P z so that the small-momentum logarithms in both C (1) L (µ 1 ) and C (1) R (µ 2 ) vanish. To extract the lightcone-DA at a specific factorization scale µ, we can relate the quasi-DA,φ(x), to the lightcone DA, ϕ(µ, x), throughφ whereĈ andV are both operators acting on the function ϕ through a convolution. Once we solve the resummed matching kernel,Ĉ RGR (µ), the lightcone DA can be extracted from the inverse matching, An easy way to implement the RGR for DA numerically is to write all the operators in the matrix representation. Both the quasi-DA and the lightcone DA are vectors on a grid of x i , ϕ i = ϕ(x i ). The matching kernelĈ and the evolution operatorV are now matrices on a grid of x i and x j ,Ĉ i j = C(x i , x j ). Thus the resummed matching in matrix representation isφ where the resummed kernelĈ RGR now is just an n × n matrix, which can be calculated row by row.
We have used an approximation of V ⊗φ = V ⊗ ϕ + O(α s ) to construct the ratio (V (0) L ⊗φ)/(V (0) ⊗φ). However, such an approximation is not applicable to all regions, because there are zero points for V (0) ⊗ ϕ. Near this point x 0 , the higher order term ∼ O(α s ) cannot be ignored, and the ratio also blows up. The point x 0 of course depends on the shape of the DA and quasi-DA curve. A test on several functions suggests that x 0 ∼ 0.2 for most DA-like functions. Thus we are still able to perform the above resummation within the region 0.25 < x < 0.75. In principle, the formalism may also work for x < 0.1 and x > 0.9.
Another issue of constructing the matching matrix is the endpoint regions of x. When either the scale µ 1 = 2xP z or µ 2 = 2(1 − x)P z becomes very small, α s (µ) increases to O(1). In this region, the perturbative expansion in α s (µ) fails, thus the matching cannot be determined through perturbative calculations. In the numerical evaluation, these rows may blow up. Since we are calculating the matching kernel in perturbation theory, the inverse matching kernel at NLO can be expressed as Applying to our matrix form, we get the inverse matching matrix which does not need the full information of the matrix to obtain a part of its inverse. So the endpoint region is no longer a problem for us to extract the resummed DA ϕ(x, µ = 2 GeV) in the intermediate x region. Those endpoint regions are thus left undetermined by perturbative calculations, and may be obtained from some nonperturbative approaches in the future. With this approach, we are able to extract the lightcone DA in the intermediate region 0.25 < x < 0.75.

Constraints from Short-distance OPE
LaMET allows us to determine the pion DA in the moderate x-region. The calculation begins to break down when x → 0 or x → 1 as shown in Eq. (3). We can, however, determine the global behavior of the pion DA from the short-distance OPE and use this information to constrain the endpoints.
At short distances, the renormalized coordinate space matrix elements can be expanded in terms of the Mellin moments of the DA where C nm (z, µ) are the Wilson coefficients, and O(z 2 Λ 2 QCD ) are higher twist effects which become relevant at distances z ≳ 0.2 fm. Fitting the matrix elements at short distances to the relevant Wilson coefficients allows us to determine the Mellin moments, ⟨ξ n ⟩. The Mellin moments describe the global information of the DA and, thus, set constraints on the endpoint region once the mid-x distribution is determined. Combined with the physical requirement that ϕ(0) = ϕ(1) = 0, it allows us to obtain the shape of the DA near the endpoints with small model dependence. This is called the "complementarity" between the large momentum expansion and the short distance factorization.
In our DA analysis, we can fit the moments from the following RG invariant and renormalization independent ratio [17]: which can be truncated at some order because higher moment contributions are negligible. Then with the Wilson coefficients known from the perturbative calculation, a fit to the short-distance ratio M(z, P 1 , P 2 ) determines the second Mellin moments ξ 2 independent of the renormalization method. Given the ϕ L (x, µ) calculated from LaMET in the mid-x region, we can model the full x-dependence, ϕ f (x, µ), as where x 0 is the minimal x we can calculate with LaMET. Then we can determine the parameter m by requiring (49) to obtain the full distribution. An alternative approach is not to constrain the endpoint region with moments, but with short distance correlations. Constructing the same full x-range distribution ϕ f (x, µ), we can first Fourier transform it to coordinate space then use the short-distance factorization in Eq. (18) to convert it to the quasi-DA correlationsh f (z, P z , µ), and fit to our renormalized matrix elementsh(z, P z , µ). This approach depends on the renormalization of our matrix elements in coordinate space, but not on the data at other momenta. We expect the second approach to give the same result as the first one. The two approaches provide a consistency check for our renormalization method with the short-distance OPE.

Numerical Results
In this work, we re-analyze the data presented in Ref. [23], measured on three lattice ensembles of lattice spacings a = {0.0582(4), 0.0888(8), 0.1207(11)} fm and pion mass m π ≈ 310 MeV generated by the MILC collaboration [51]. The analysis starts with the same bare matrix elements extracted from a two-state fit to the lattice correlators.

Renormalization
The method used in Ref. [23] for renormalization was the regularization-independent momentum subtraction (RI/MOM) scheme. However, this method has some problems when dealing with the linear divergence in the nonlocal operator with a spatial Wilson line, as well as generating unknown nonperturbative effects at large distances [52]. We deal with these issues by working in the hybrid scheme [29] with the LRR-improved ratio scheme at short distances from the self-renormalization [30], as discussed in the previous section. The renormalization factors Z R (z, a) at short distances are obtained from P z = 0 matrix elements of the pion PDF to remove the linear divergence, where ∆I is a conversion constant in different schemes, adjusted through fitting to make sure the small z-correlations are matched to MS results. The fitting parameter d = −0.53 represents the NLO RG evolutions on the lattice. The term k a ln(aΛ QCD) is the linear divergence with fitting parameters k and Λ QCD which are not uniquely determined due to the intrinsic ambiguity. By choosing a set of fitting parameters k = 3.3 fm −1 GeV −1 and Λ QCD = 0.1 GeV [30] as scheme τ ′ , we can determine a corresponding m 0 (τ) = 0.151 GeV [32] to relate the P z = 0 lattice matrix elements to the perturbative calculation of C LRR 00 (z, µ, τ), as defined in Sec. 2 with the asymptotic form, which eliminates such an ambiguity. The term f (z)a as a fit parameter incorporates the discretization effects and the remaining terms come from the resummed logarithmically divergent dependence on a.
After removing the linear and logarithmic UV divergences through Eq. (4), we are able to extrapolate to the continuum limit to take out additional discretization effects at finite P z , which is a simple process of fitting the renormalized matrix elements at different lattice spacings but at a fixed z value to a linear function: for some function c(z) whereH R is defined in Eq. (22). We carry out this extrapolation for a continuous curve after interpolating our data on three lattice spacings. The matrix element after renormalization and continuum extrapolation  is shown in Fig. 1, which shows a good consistency among different lattice spacings. A comparison of renormalized matrix elementsH R at different momenta with the LRR perturbative result at P z = 0 is shown in Fig. 2. We can clearly see that the distribution approaches the P z = 0 perturbative calculations C LRR 00 (z, z −1 ) when the momentum decreases. Then a ratio to the P z = 0 perturbative results for DA is taken at short distance |z| < z s to convert to the hybrid scheme, whereH R is defined in Eq. (22). The P z > 0 matrix elementsH R (z, P z ) are obtained from our data renormalized by Eq. 51, while theH(z, P z = 0) ∝ P z = 0 vanishes in DA measurement, thus we use the Wilson coefficient C 00 from the perturbative calculation in MS scheme. The parameter z s that separates the nonperturbative and perturbative regions must be much larger than the lattice spacing to avoid discretization effects but not so large as to necessitate higher-twist terms in the OPE. In our calculations, we choose z s = 0.18 fm. Note that there is no modification of z dependence in the second term to avoid introducing unwanted non-perturbative effect. Figure 3 shows a comparison between the ratio M hybrid (z, P z , 0) obtained from fixed-order self-renormalization and the LRR-improved self-renormalization. An obvious problem in the fixed-order renormalization is the hump near z = 0.1 fm, which suggests a negative second moment ξ 2 < 0, irreconcilable with the OPE predictions at short distance. The LRR-improved renormalized matrix elements, on the other hand, show good consistency with OPE predictions at short distances. This comparison demonstrates that the modification from LRR is necessary for a correct renormalization. Both the continuum extrapolation and the conversion to the hybrid scheme are linear, thus the two steps commute with each other.

Extracting x-dependence of quasi-DA
In order to Fourier transform our coordinate space correlations to momentum space and extract the x-dependence, we need to extrapolate our matrix elements to infinite distance. We first convert our position-space variable to quasi light-cone distance λ ≡ zP z as introduced in Sec. 3.1. Since P z is fixed for a single calculation, large distance corresponds to large λ. Although the large-λ correlation becomes extremely noisy at a finite momentum from the lattice, which in principle makes it impossible for us to know the longtail information, the distribution is not arbitrary. A general consideraction of coordinate space correlations suggests an algebraic decay along with an exponential decay [29,25]. The constraints allow us to reduce the error in the large-λ region and extract the x-dependence of the quasi-DA.
Based on these constraints, we can extrapolate our matrix elements in position space to the corresponding inverse Fourier transform [29,25]:H where λ 0 is a large constant describing the correlation length and depends on the hadron momentum, and the terms (c 1 , n) are fitting parameters. Note that at long distances in our hybrid scheme, the ratio M hybrid (λ, P z ) only differs fromH R (λ, P z ) by a constant factorH R (z s , 0), so they have the same functional form. We then fit the longtail of M hybrid (λ, P z ) to Eq. (54). An example of the large-λ extrapolation is shown in Fig. 4.  With a full-range coordinate space correlation, we are able to extractφ(x, P z ) through a Fourier transformation in our chosen renormalization scheme,φ as shown in Fig. 5. Although we used a model to describe the large-λ behavior, we should address that the final result is not sensitive to the highly suppressed long tails. To illustrate this, we use two different models corresponding to Eq. (54), one is not to include the exponential decaying factor (labelled as "Power Decay"), the other is by fixing the correlation length λ 0 = 50 (labelled as "Exp Decay2") to check the sensitivity of DA to the long-tail model assumptions. The comparison is shown in Fig. 6, showing small discrepancies in mid-x region when compared to the statistical error in blue band. The endpoint regions are more sensitive to the long-tail modeling, but we only calculate the mid-x region of light-cone DA directly and the endpoint regions are obtained from modeling. So the long-tail modeling dependence has little influence on our final determination of the full-x distribution.

Matching to obtain DA in mid-x region
After obtaining the quasi-DA in momentum space, we can then apply the matching to obtain the lightcone DA. Firstly, we apply the fixed-order matching kernel, modified with LRR, at scale µ = 2 GeV without the large log resummation. The fixed order matching appears to be valid for the full x region. However, as we discussed, it is just an artificial effect. As we approach the endpoints, the higher-order large logs can no longer be neglected and have to be resummed. We thus apply the matching with RGR, and show the comparison in Fig. 7. As we discussed in    We also find that the result is insensitive to which LRR method is used, and the scale choice for RGR by changing the initial scale of RGR from 2xP z to 2cxP z with c ∈ [0.75, 1.5], suggesting only < 3% difference, as shown in Fig. 8.  They are all consistent within error, up to < 3% difference.

Full x-dependence for lightcone DA
Now we have the distribution determined for the mid-x region, while the endpoint regions are still unknown from the LaMET approach. Fortunately, in coordinate space, the matching coefficients in small-z region is perturbative, thus can be resummed safely. Applying these matching coefficients to lattice data, we are able to obtain the lightcone correlation in a certain range of correlation length λ = zP z , which contains the global information of the x-dependent DA, such as its moments. With the information from the mid-x region, we can complete our picture of the extracted DA by utilizing the small-z information to constrain the endpoint behavior, as suggested in Ref. [34].
Near the endpoints, we can parametrize the DA as a power of x or 1 − x as in Eq. (48). To ensure continuity, we require that the parametrized form coincide with our mid-x results at x = x 0 . We convert this parametrized DA into coordinate space, apply the short-distance matching, and fit the result to our renormalized matrix elements. Figure 9 shows the comparison from the parametrized DA and our lattice data at short distances, which suggests good consistency. Besides that, we allow some model-dependence by adding a small correction, where sin(b) is to guarantee that the size of this correction term is not so large as to cause a sharp turn at the junction point, i.e., the different regions are smoothly connected. The same modification is symmetrically applied to x → 1 − x, and a model-dependence is included as the systematic error by calculating the difference between the modified model and the original one in Eq. (48). The final estimation taking into account such a systematic error is shown in Fig. 10. We can estimate the moment from the full x-dependence from Eq. (45) with Eq. (48). We get ⟨1⟩ =0.999(5), ξ 2 = 0.306 (19), (57) which are in good agreement with the theoretical normalization ⟨1⟩ = 1, and the second moment ξ 2 = 0.298 (39) obtained from the renormalization-independent OPE fit to Eq. (47). This self-consistency is a strong support for our renormalization method. In Fig. 11 we compare with previous model-dependent calculations and lattice results, including the Dyson-Schwinger Equation (DSE'13) [53], the prediction of the light-front constituent-quark model (LFCQM '15) [54], the OPE reconstruction from local second moment calculations (RQCD'19) [12], the lattice calculation from LPC (LPC'22) [25], and the reconstruction from fitted moments by ANL/BNL collaboration (ANL/BNL'22) [17]. Our final result suggests the flattest and broadest distribution among all these calculations, as we can also tell from the large second moment in our data. This may have a big impact on the phenomenology of pion hard exclusive processes. At large Q 2 , the πγγ * transition form factor F πγγ (Q 2 , 0) [4,3,5,55] and the pion electromagnetic form factor F π (Q 2 ) [56] are both sensitive to the shape of the DA ϕ(x, Q 2 ). In general, since both are enhanced near the threshold x → 1 [57,58], a broader DA will predict the form factors to be larger at large Q 2 . However, the factorization of the exclusive processes are known to be problematic near the threshold [47,59,60], and the scale setting in the pion electromagnetic form factor F π (Q 2 ) also causes large uncertainty to its estimation from the DA [61,57]. Due to these complications, the study of these phenomenologies are beyond the scope of this work, but a more detailed and systematic study is needed in the future to completely understand the impacts.

Conclusion
In this paper we have computed the pion distribution amplitude with momentum fraction in the range x ∈ [0, 1] with improved handling of three sources of systematic errors: removing the O(Λ QCD /xP z ) power correction from intrinsic ambiguities, resumming the small-momentum logarithms, and constraining the distribution near the endpoints from short distance correlations. We renormalize the matrix elements in the hybrid scheme with the LRR-improved self-renormalization factors at short distance. Then an LRR-improved matching kernel is used, along with a two-scale resummation, to obtain the pion DA in the mid-x region x ∈ [0.25, 0.75]. We then model the endpoint region with a power law function, allowing a small variation, to reconstruct coordinate space correlations and fit to our data. The second Mellin moment determined from the full-x dependence was ξ 2 = 0.302 (23) and from the renormalizationindependent short-distance OPE was ξ 2 = 0.298 (39). These two results are in good agreement and give us confidence in the determination of the endpoint region of the DA. Our final result suggests a broad distribution of the pion DA. It has the potential for a big impact on the form factors of the DA at large Q 2 , and will be investigated in detail in the future.