Lattice QCD calculation of the pion distribution amplitude with domain wall fermions at physical pion mass

: We present a direct lattice QCD calculation of the x -dependence of the pion distribution amplitude (DA), which is performed using the quasi-DA in large momentum effective theory on a domain-wall fermion ensemble at physical quark masses and spacing a ≈ 0 . 084 fm. The bare quais-DA matrix elements are renormalized in the hybrid scheme and matched to MS with a subtraction of the leading renormalon in the Wilson-line mass. For the first time, we include threshold resummation in the perturbative matching onto the light-cone DA, which resums the large logarithms in the soft gluon limit at next-to-next-to-leading log. The resummed results show controlled scale-variation uncertainty within the range of momentum fraction x ∈ [0 . 25 , 0 . 75] at the largest pion momentum P z ≈ 1 . 85 GeV. In addition, we apply the same analysis to quasi-DAs from a highly-improved-staggered-quark ensemble at physical pion mass and a = 0 . 076 fm. By comparison we find with 2 σ confidence level that the DA obtained from chiral fermions is flatter and lower near x = 0 . 5 .


Introduction
Understanding the structure of pions is of special importance as they are the lightest pseudo-Nambu-Goldstone bosons of chiral symmetry breaking in strong interactions.One particular quantity of interest is the pion distribution amplitude (DA) ϕ(x) which describes the probability amplitude of finding pion on the light-cone in a quark-antiquark pair Fock state, each carrying a momentum fraction of x and 1−x.The rich phenomenology of the pion DA originates from its universality as inputs to exclusive processes and form factors at large momentum transfer Q 2 ≫ Λ 2 QCD .For example, the scattering amplitudes of semileptonic B-meson decay [1,2] and deeply virtual meson production processes [3], which have been used to probe physics beyond the Standard Model and to extract the generalized parton distributions, respectively, are proportional to convolutions involving the pion DA in the x-space.Therefore, knowing the x-dependence of DAs is key to making predictions for the hard exclusive processes.However, they are only weakly constrained by experiments so far [4][5][6][7], which makes it highly desirable to calculate their x-dependence from first-principle methods like lattice QCD.
The pion DA ϕ(x, µ) is defined from a light-cone correlation as where f π is the pion decay constant, W (0, η − ) = P exp −ig η − 0 ds n µ A µ (ns) is the Wilson line between the two light-cone coordinates 0 and η − = (η 0 + η 3 )/ √ 2 with P denoting the path-ordering operator.Due to the real-time dependence of light-cone correlations, it is impossible to directly calculate the full pion DA on a Euclidean lattice.
There are different approaches to extract information of the pion DA from lattice QCD, including the calculation of its moments ⟨(2x − 1) n ⟩ [8][9][10][11][12][13], short distance factorization (SDF) of nonlocal correlations [14][15][16][17][18][19], and the large momentum effective theory (LaMET) [20][21][22][23][24][25][26][27].The calculation of moments is based on the light-cone operator product expansion (OPE) of the pion-DA correlator, which encodes the DA moments in twist-two local operators [8].The bottleneck of this method is the increasing noise in higher moments and power-divergent operator mixings for n > 3.In the SDF approach, the spatial correlation function in the pion state can have an OPE [28] or be factorized into the convolution of the light-cone correlation and a perturbative matching kernel at a distance |z| ≪ Λ −1 QCD .Therefore, this method can provide information of the first few moments, which in principle can go beyond n = 3, or the light-cone correlation within a limited range.However, it still requires a model assumption of the shape of DA to reconstruct its x-dependence [16,19].The LaMET approach provides a direct calculation of the x-dependence through an effective theory expansion in the momentum space, up to power corrections suppressed by the parton momenta xP z and (1 − x)P z , which allows for a reliable calculation of the DA within a moderate range of x.
Since the first LaMET calculation of the pion DA [23], various lattice artifacts and theoretical systematics have been studied.One of the key systematics is the lattice renormalization of the quasi-DA correlators, which suffer from the linear power divergence that must be subtracted at all distances [29][30][31][32].The regularization-independent momentum subtraction (RI/MOM) [33][34][35][36] and ratio [37,38] schemes were proposed to cancel the linear and other ultraviolet (UV) divergences, but both introduce uncontrolled non-perturbative effects at large distance for LaMET calculation.Then, the hybrid scheme [39] was proposed to overcome this problem with a subtraction of Wilson line mass at long distances [40,41], but there is still a remaining linear renormalon ambiguity of order Λ QCD that contributes to a linear power correction in the LaMET expansion.Recently, this issue was resolved with the leading-renormalon resummation (LRR) method [27,42], which removes such an ambiguity in lattice renormalization and LaMET matching, thus improving the power accuracy to sub-leading order.Apart from lattice renormalization, there has also been significant progress in the perturbative matching, from the next-to-leading order (NLO) matching kernel for quasi-DA [23,24,[43][44][45] to the renormalization group resummation (RGR) [46,47] and threshold resummation [46,48,49].On the numerical side, the first continuum extrapolation of the DA was carried out in Ref. [25] with the RI/MOM scheme at unphysical pion masses, which observed a sensitive dependence on the pion mass.This calculation was succeeded by a continuum extrapolation at physical pion mass with NLO hybrid-scheme renormalization and matching, but without LRR [26].The LRR method was first implemented in Ref. [27] with NLO matching and RGR.Meanwhile, the pion DA moments have also been calculated using the SDF approach at NLO accuracy [16,19].Nevertheless, no lattice calculation in the literature has included the threshold resummation yet.
So far, all the lattice calculations of pion DA x-dependence have been performed with highly-improved-staggered-quark (HISQ) fermion actions [50,51] or clover fermion actions [52], both of which explicitly break the chiral symmetry.Since the pions are the pseudo-Nambu-Goldstone bosons of chiral symmetry breaking in QCD, it is then very interesting to study how the chiral symmetry of lattice action affects their structure such as the DA.The overlap fermion [53][54][55][56] and domain wall fermion (DWF) actions [57][58][59] are known to preserve chiral symmetry on the lattice, although these ensembles are more expensive to generate.Therefore, performing the same calculation on chiral fermion ensembles will provide us with knowledge about the effects of chiral symmetry breaking on the pion structure.
In this work, we present the first direct calculation of the x-dependence of pion DA on a lattice ensemble with the DWF fermion action at physical pion mass and spacing a ≈ 0.084 fm.We renormalize the bare pion quasi-DA matrix elements with the hybrid scheme, and then match them to the continuum MS scheme with LRR.After Fourier transforming to the x-space, we match the quasi-DA to the light-cone with next-to-next-to-leading logarithmic (NNLL) threshold resummation and NLO matching.For the first time, we reformulate the recently developed threshold resummation technique for the quasi-PDF case [48] to work for the quasi-DA, which improves the estimate of systematic uncertainties near the end-point regions x → 0 and x → 1.With the reformulated threshold factorization, we also find out that the complicated two-scale resummation [27] of the Efremov-Radyushkin-Brodsky-Lepage (ERBL) logarithms [60][61][62][63] is now factorized into two separate pieces with each involving only a single physical scale, thus making it straightforward to understand and implement the RGR.Finally, we observe a notably flat DA in the region x ∈ [0.25, 0.75], while the results beyond this range become unreliable due to the breakdown of perturbation theory.Moreover, compared with the same analysis applied to the HISQ data [19] at physical pion mass and lattice spacing a = 0.076 fm, we observe that the pion DA from the DWF ensemble is slightly flatter near x = 0.5 than that from the HISQ ensemble.Our results suggest that the pion structure is not sensitive to chiral symmetry on the lattice.This work is organized as follows.In Sec. 2, we present our lattice setup of the calculation, and show how the raw lattice data are processed to extract the matrix elements of pion DA.In Sec. 3, we analyze the bare matrix elements to extract the DA moments using the SDF or OPE approach, and discuss how the matrix elements are properly and consistently renormalized in the hybrid scheme with LRR to obtain x-dependent quasi-DA.In Sec. 4, we derive the formalism to resum both ERBL and threshold logarithms in the perturbative matching and implement the improved matching to extract the light-cone DA, and compare with the same analysis applied to data on a HISQ ensemble [19].Finally, we conclude in Sec. 5.

Lattice set up
In this calculation, we used a 2+1-flavor domain-wall gauge ensemble generated by RBC and UKQCD Collaborations of size N 3 s ×N t ×N 5 = 64 3 ×128×12, denoted by 64I [64].The quark masses are at the physical point and the lattice spacing is a = 0.0836 fm.55 gauge configurations were used in this calculation.The quark propagators are evaluated from Coulomb-gauge-fixed configurations using deflation based solver with 2000 eigen vectors.
In addition, we used the boosted Gaussian momentum smearing [65] to improve the signal.The Gaussian radius was set to be r G = 0.58 fm.We chose the quark boost parameter j z to be 0 and 6 [66,67] which are optimal to hadron momentum P z = 2πn z /(N s a) with n z = 0 and 8, so that the largest momentum in our calculation is P z = 1.85 GeV.Since only two-point functions are involved in this calculation [19], measurements at other momenta (n z ∈ [0, 3] for j z = 0 and n z ∈ [4,8] for j z = 6) were also computed through contractions using the same profiled quark propagator.To increase the statistics, we ultilzed the All Mode Averaging (AMA) technique [68] with 2 exact and 128 sloppy sources for momenta n z ∈ [4,8], while 1 exact and 32 sloppy sources for n z = [0, 3].The tolerance of the exact and sloppy sources are 10 −8 and 10 −4 respectively.To suppress the ultraviolet (UV) fluctuations and enhance the signal-to-noise ratio of the matrix elements with long Wilson links, we employed Wilson flow [69], with a flow time t F = 1.0 (roughly corresponds to a smearing radius √ 8a 2 ).Utilizing the symmetry of the data, we further average the forward and backward correlators at Euclidean time slice τ and N t −τ , as well as averaging the quark bilinear separation z and −z, with the corresponding parity.In total, we have effectively 28,160 measurements for the largest momentum P z = 1.85 GeV.
In accordance to the quasi-DA definition we measure three different types of correlators, C ππ , C πO 0 , and C πO 3 , to extract the quasi-DA matrix elements: where O π = ψγ 5 ψ is the smeared pseudo-scalar interpolator, that has an overlap with pion ground state is the spatial Wilson line that keeps the non-local operator gauge invariant.The Gaussian momentum smearing is applied to the pion interpolator O π (x, t), so C ππ is smeared at both the source and the sink, while C πO 0 and C πO 3 are smeared only at the source.Note that in general there could be a mixing between the non-local operators ψγ 5 γ t W z ψ and ψγ x γ y W z ψ [33,70], but such a mixing is proportional to the explicit chiral symmetry breaking, thus vanishes specifically on the DWF ensemble.
By expanding these correlators in a tower of energy eigenstates, we can extract the energy spectrum and the coefficients, with the ground-state coefficients where H γ t/z γ 5 (z) are the matrix elements of pion quasi-DA, which is normalized to H γ t/z γ 5 (0) = 1.The P z = 0 correlators are non-vanishing for C ππ (t) and C πO 0 (t, z), so they could be used to extract the pion mass.
The effective mass for the C ππ correlator at different momenta are shown in Fig. 1.At small Euclidean time, the excited states effect is important, thus making the effective masses of the smeared-smeared C ππ correlator larger than the actual ground state energies.At large Euclidean time, the effective masses decays to approach ground state plateaus, which are consistent with the dispersion relation E(P z ) = P 2 z + m 2 from zero momentum correlators, plotted as colored lines.Thus we two-state fit the energy spectrum at z = 0 using the dispersion relation as a prior for E 0 .The first-excited state of the pion using smeared correlators on the lattice has been studied in previous works [66,71].In these works it was suggested that the first excited state is the π(1300) state.Therefore, the energy of the first excited pion state for different moment can be estimated using the dispersion relation and the mass of π(1300) state.We use this estimate as a prior for E 1 , with a width of 0.5 GeV.The first-excited state energies from two-state fit are shown in Fig. 2 as function of t min , with t min being the lower limit of the fit interval.By including the excited state contribution, we are able to utilize data at smaller Euclidean time, where the plateau of effective mass has not been reached, and still get good fit quality.Two examples of the fits at n z = 8 and z = 0 and z = 3a with different fit ranges are shown in Fig. 3. To compensate the exponential decay in Euclidean time, we multiply the correlators with a factor of e Ē0 t .There is a good consistency between data and our t min = 4a fit bands.The fitted ground state E 0 and the first excited state E 1 are shown in Fig. 4 as a function of n z .These results are consistent with the priors based on the dispersion relation, depicted as blue and orange bands.Once E 0 and the ground state coefficients A O 0 0 and A π 0 in Eq. (2.3) are known, we can calculate bare f π , Although we have not calculated the renormalization constant Z A and Z S to obtain the renormalized f π , which makes it infeasible to directly compare with the physical value, the consistency of f bare π from different fits can be used as a criteria of the fit quality besides the χ 2 /d.o.f and p-value.The f bare π calculated from the coefficients of these different fits are shown in Fig. 5.The plot suggests that the fitted f π is consistent for different momenta with different choices of t min , although the error increases significantly for large pion momenta.Taking energies from the above two-state fit of local correlators as inputs, we fit the non-local correlations |z| > 0 to extract the coefficients A 0 (z) in Eq. (2.4).The bare matrix elements of C πO 0 and C πO 3 are both shown in Fig. 6.Note that the energy spectrum are fitted separately, in order to obtain a better fit quality in each data set.In our fit results, both imaginary parts are consistent with zero, because the quark and anti-quark are symmetric in the pion.Interestingly, although we can see from Eq. (2.4) that the coefficient A O 3 0 (z) should statistically vanish at P z = 0, after normalizing the matrix elements as 0) sample by sample, the ratio averages to non-zero values.So we are able get P z = 0 matrix elements for C πO 3 although with large errors.
3 Extraction of pion quasi-DA

Moments of pion DA
The quasi-DA correlator has been proven to be multiplicatively renormalizable [30][31][32],  where the renormalization factor Z(z, a) includes the linear and logarithmic UV divergences.Since the UV divergences are exactly the same for the same z at any P z , we can cancel Z(z, a) by taking the ratio of matrix elements at two different momenta [37,38], The above ratio is renormalization group invariant, so we suppress the µ dependence in the MS renormalized correlation function H R (z, P z ), which can be factorized at short-distance which depends on the non-perturbative moments of the light-cone-DA, and the perturbative matching coefficients C nm (z, µ).It allows us to extract the moments ⟨ξ m ⟩(µ) from the ratio M(z, P 1 , P 2 ) when the higher twist correction are suppressed For the symmetric pion DA, only even moments of ξ are non-vanishing.We can expand Eq. (3.5) to m, n = {0, 2, 4} at NLO [19,72]: where i = 0 or 3, L = ln z 2 µ 2 e 2γ E

4
. The triangular matrix C indicates a non-multiplicative renormalization group evolution [63].The operator itself follows the renormalization group (RG) equation where γ C is the anomalous dimension of the quark bilinear operator that has been calculated up to 3-loop order [73].Besides, note that which can be derived from the fact that only the n = m = 0 term remains in the P z = 0 matrix element.Moreover, Eq. (3.7) must be satisfied for each individual term in n, so there is certain cancellation of the µ-dependence between the coefficients and moments.The 1-loop evolution can be directly read from the log terms of C γ i 00 − C γ i nm in Eq. (3.6): Beyond 1-loop order, the ERBL kernel has been calculated in momentum space [74], and γ nn is the same as those of PDF operators.The off-diagonal part of γ nm has been calculated up to 3-loop order using conformal symmetry [75].Here we quote their number and convert them into the Mellin basis, for in units of α 2 s .So we can either fit the moments {⟨ξ 2 (µ)⟩, ⟨ξ 4 (µ)⟩} with NLO Wilson coefficients, or fit with RG-resummed (RGR) Wilson coefficients.In the latter case, we fit the moments at an initial scale ⟨ξ 2 (µ i = 2e −γ E z −1 )⟩, where the log terms in C nm (z, µ i ) vanish, then evolve to µ = 2 GeV by solving Eq. (3.9).We examine the scale variation by choosing 2} to estimate the uncertainties from higher-order perturbation theory.The fitted ratio and the extracted moments are shown in Fig. 7.In the left figure, the difference between the fitted ratios from the NLO and NLO+RGR Wilson coefficients are negligible, because the two different fits are eventually optimized to the same function of (zP z ) 2 that best describes the data.The moments are then solved from these same coefficients of (zP z ) 2n with different Wilson coefficients, thus they could be quite different.We show the fitted moments with both statistical and scale variation error bars.Note that the scale variation becomes large when z increases, because the scale z −1 / √ 2 ∼ Λ QCD becomes non-perturbative.The fact that the fit results depend on the z-value of the data indicates a non-trivial discretization effect, and it could only be determined accurately when multiple lattice spacings are included, or on fine lattice spacings where we have more reliable data points z ≫ a while still stays in the perturbative region z ≪ Λ −1 QCD .

Renormalization
In the previous section, we have extracted the moments of pion DA through a renormalizationindependent approach.To directly calculate the x-dependence ϕ(x, µ), however, we have to first properly renormalize the correlation functions.Among multiple renormalization schemes, the hybrid renormalization scheme [39] is preferred in quasi-DA calculation because it allows a perturbative matching to MS scheme at all z, where the factorization theorem of quasi-DA is proven.In the hybrid scheme renormalization, where δm(a) ∼ a −1 is the mass renormalization parameter that removes the linear divergence, and the short-distance renormalization factor Z R (z, a) can be some lattice matrix elements with the same divergence, for example, the same operator measured in other external states.A widely-used choice is the matrix element in the same hadron state at rest [37], At long range |z| > z s , the hybrid scheme requires the determination of the mass renormalization parameter δm(a).However, its determination contains an infrared ambiguity of ∼ O(Λ QCD ) [39], which can be labelled as a regularization scheme dependence δm(a, τ ).When combined with the renormalon ambiguity in the perturbative matching kernel C(x, y, µ, P z ), it results in a linear correction in 1/xP z in the final extraction of light-cone DA [39].It is proposed to absorb the linear ambiguity into a non-perturbative parameter m 0 independent of the hadron momentum [39,40,42], such that m 0 can be extracted from P z = 0 matrix elements, and is applied to the renormalization of P z > 0 data by replacing δm → δm(a) + m 0 in Eq. (3.11).The value of m 0 is still ambiguous and z-dependent, except when both ambiguities in the extraction of δm(a, τ ) and the perturbative matching kernel C(x, y, µ, P z , τ ) are regularized in some specific scheme τ .After removing the linear divergence in a specific regularization scheme, the z-dependence of short distance z ≪ Λ −1 QCD matrix elements should be well described by perturbation theory, up to power corrections and lattice artifacts.Since we have only one lattice spacing, we ignore the z-dependent discretization, and extract m 0 (τ ) by requiring that the P z = 0 renormalized matrix element matches the perturbatively calculated Wilson Coefficient C 0 (z, µ, τ ) with renormalon regularized in scheme τ , where is the renormalization group evolution factor that cancels the renormalization scheme dependence between the lattice scheme H B (z, P z = 0, a) and the MS scheme result C 0 (z, µ, τ ) [42], in which the anomalous dimension γ(α) have been calculated to 3-loop order [73] and the beta function β(α) have been calculated to 5-loop order [76,77] in MS scheme.Here we apply the leading renormalon resummation (LRR) method as introduced in Ref. [42] to regularize the linear ambiguity and extract m 0 (τ ), such that the linear power corrections are cancelled when a corresponding LRR-improved matching C LRR (x, y, µ, P z , τ ) is applied.Note that we have only one single lattice spacing a = 0.0836 fm, so with a and τ fixed, δm(a, τ ) is just a constant.Thus we redefine δm ≡ δm(a, τ ) + m 0 (τ ) and extract it together from P z = 0 data via Eq.(3.13).Taking the ratio between two adjacent z's, we get The regularization scheme τ of the Wilson coefficient C 0 (z, µ, τ ) ≡ C LRR 0 (z, µ) is defined as an LRR improved perturbation series with principal value (PV) prescription.At NLO, its form is [42] for the operator O i , where b = β 1 /2β 2 0 and c 1 = (β 2 1 −β 0 β 2 )/(4bβ 4 0 ) are from higher orders in the QCD beta function, N m (n f = 3) = 0.575 is the overall strength of the linear renormalon estimated from the quark pole mass correction [78,79].To stay in the perturbative region and avoid higher twist contribution, in principle we should not go beyond z = 0.3 fm.Also note that z = a usually suffers from non-trivial discretization effects, as we will show in a later section.Thus we use z = {2a, 3a} to extract δm.For comparison, we also used fixed-order result C FO 0 (z, µ) = C LRR 0 (z, µ) Nm=0 , and the RGR result ) in Eq. 3.14 to extract δm.The scale variation is examined by vary the initial scale of RGR µ = 2cz −1 e −γ E with c = {1/ √ 2, 1, √ 2}, and we show the comparison for different extractions in Fig. 8.The LRR method significantly improves the scale variation and convergence of perturbation theory.Also, the δm obtained from two different operators are consistent, suggesting that our extracted δm is universal for this Euclidean non-local quark-bilinear operator in a pion external state.Since the normalized matrix elements H B γzγ 5 (z, P z = 0) are very noisy, we take the result δm = 0.397 +0.032 −0.008 GeV extracted from H B γtγ 5 (z, P z = 0) at NLO.Then the same δm is used to renormalized both H B γzγ 5 and H B γtγ 5 .After removing the linear divergence and regulating the renormalon ambiguity, the matrix elements can be compared with OPE reconstruction if the first few moments are given.Taking the moment we fitted in Sec.3.1, we show the comparison of the matrix elements e z(δm) H B (z, P z , O) and RG-invariant ratio M (z, P 0 , P z , O) with their OPE reconstruction Eq. (3.3) and Eq.(3.5) in Fig. 9, with RGR and LRR corrections to the Wilson coefficients.We also tried to introduce higher moments in the OPE reconstruction, and found their contribution to be sub-percent level in this regime.Here we have used P 0 = 0 for O 0 and P 0 = 0.23 GeV for O 3 as the denominator when constructing the ratio.There is a good agreement for the RG-invariant ratio, but a noticeable overall deviation for renormalized matrix elements e z(δm) H B (z, P z , O 0 ) at z = a.This large deviation indicates non-negligible discretization effects at z = a, which appears to be universal among all momenta.When taking a ratio, such discretization effects are cancelled between two different momenta, thus the ratio is more consistent with OPE reconstruction.Nevertheless, at z = {2a, 3a}, the matrix elments are still roughly consistent with OPE reconstruction, suggesting that the renormalization is done correctly.The consistency for both operators also suggests that δm is not sensitive to the Dirac structure in the operator.Note that the discretization effects in z = a matrix elements for O 3 are not as large as O 0 .So if we take a ratio of different operators, extra discretization effect will be introduced by H B (a, 0, γ 5 γ t ).Thus when implementing hybrid scheme renormalization, H B (z, 0, γ 5 γ t ) is only used to renormalize O 0 .Considering the good consistency between H R (z, P z , γ 5 γ z ) and OPE reconstruction, we assume that the discretization effect at z = a is much smaller for γ 5 γ z , thus we can use perturbatively calculated H B (z, 0, γ 5 γ z ) to renormalize it, according to Eq. (3.13), where µ 0 = 2e γ E z −1 , and I lat (a) is the evolution factor in lattice scheme, and is a constant for fixed lattice spacing that can be tuned to make sure Z R (z < z s ) is consistent with H B (z, 0, γ 5 γ z ).The renormalized matrix elements in the hybrid scheme Eq. (3.11) are shown in Fig. 10.We find that the matrix elements at different momenta saturate to a universal shape at large zP z , except for scaling violation in z 2 or 1/P 2 z which is not distinguishable from the statistical uncertainties.The non-smoothness of the renormalized matrix elements near z s in Fig. 10 is the nature of hybrid scheme, and has been taken into account in the perturbative matching kernel.We discuss this feature in more details in Appendix A. which is not distinguishable from the statistical uncertainties.

x-dependent quasi-DA
With the LaMET approach, we aim at calculating the local x-dependence directly in a certain range with controlled systematics.On the other hand, lattice calculations of the matrix elements are naturally performed in the coordinate space.Therefore, a Fourier transformation is needed to obtain the x-dependence of quasi-DA first.Due to the worsening signal-to-noise ratio as z increases, lattice data are limited to a certain range of z or zP z , which makes an exact Fourier transform impossible.However, although the long tail of the correlation is unknown, it follows certain physical constraints that prevent it from going out of control [39].By performing an extrapolation in zP z using such constraints, one can reduce the uncertainties in the long-tail region, thus eliminating unphysical oscillatory behaviors in the x-space after the Fourier transform.One strong physical constraint on the long-tail distribution is the finite correlation length λ 0 ∝ P z in the coordinate space, which results in an exponential decay exp{[−λ/λ 0 ]} of the Euclidean correlation functions at large λ = zP z .For example, we can model the long-tail region to the following form [39], where the parameterization inside the round brackets is motivated from the Regge behavior [80] of the light-cone distribution near the endpoint regions.For symmetric pion DA, To estimate the model dependence, we perform three different fits of the long tail, considering where the correlation functions start to be dominated by the exponentially decay, or more conservatively, assuming no exponential decay at all.These attempts include • fitting data from λ ≈ 8.5 to Eq. (3.18), labeled as "Exponential"; • fitting data from λ ≈ 8.5 to Eq. (3.18) without the exponential decay e −λ/λ 0 , labeled as "Algebraic"; • fitting data from λ ≈ 11 to Eq. (3.18), labeled as "Large λ".
The results are shown in Fig. 11.The long-tail extrapolation will eventually introduce about ±5% systematic uncertainties to the light-cone DA near x = 0.5 in our case.To reduce this uncertainty, a more precise measurement of the long-range correlation of quasi-DA is necessary, especially for determining the oscillating behavior of the long tail.4 Matching to the light-cone DA The light-cone DA is extracted from quasi-DA through an inverse matching with power corrections starting from quadratic order in Λ 2 QCD /P 2 z after the linear corrections are eliminated by the LRR improvement: with the perturbative matching kernel C(x, y, µ, P z ) calculated up to NLO [45], and in the hybrid scheme [39], where x = 1 − x, and the plus function defined in a certain range [a, b] is For γ z γ 5 operator, there is an additional correction term C γzγ 5 = C γtγ 5 + ∆C γzγ 5 , (4.5) Two important systematics have been discussed in Ref. [27], including the leading power correction and the large small-momentum logarithmic contributions.The linear power correction has been demonstrated to be eliminated through LRR [42], so we apply the same correction on the matching kernel here.On the other hand, the small-momentum logarithmic contributions can be resummed [47] by solving the (ERBL) RG equation, but there is no well-established method due to the more complicated two-scale nature of the DA matching kernel [27].
Note that the logarithm related to quark (antiquark) momentum fraction in Eq. ( 4.2) has the form x y ln x or x ȳ ln x.The coefficients of the logarithm are thus suppressed when x → {0, 1}, whose contributions vanish when the quark (antiquark) momentum becomes very soft, except for the limit when x → y.Therefore, the ERBL logarithm only needs to be resummed in the threshold limit x → y.And as we will show below, the resummation of small-momentum logarithms will also become more straightforward in the threshold limit.

Physical scales in the matching kernel: threshold limit
It has been shown that there are two physical scales in the quasi-DA matching, corresponding to the quark momentum 2xP z and antiquark momentum 2(1 − x)P z , respectively [27].These two scales becomes apparent when considering the threshold limit x → y of C(x, y, µ, P z ), where it can be factorized into a heavy-light Sudakov form factor H(xP z , xP z , µ) and a jet function J(|x − y|P z , µ) [48,81], up to higher orders of (y − x) (note that the leading term is singular in (y − x), thus the corrections starts from the 0-th order), In coordinate space, the factorization is multiplicative.Thus the order of the Sudakov factor H and the jet function J in the momentum space factorization Eq. (4.6) does not matter in the threshold limit.The heavy-light Sudakov factor is obtained by matching the heavy-light currents to soft-collinear effective theory [82], which in our case is the product of two components from two heavy-light currents in coordinate space [82,83], where the signs depend on whether the momentum of the quark (antiquark) is aligned with (opposed to) the wilson line direction.For example, From now on we will focus on the physical region 0 < x < 1 only, because the non-physical region is always far away from the threshold limit x → y when 0 < y < 1.The heavy-light Sudakov factor in the quasi-DA differs from the quasi-PDF case in two aspects: the two heavy-light vertices have different momenta in quasi-DA; their phases are opposite in the quasi-DA.Although the Sudakov factor in Eq. (4.7) explicitly depends on the momentum fraction x, its imaginary part has z-dependence through sign(z) of the Wilson line direction.When Fourier transformed to the momentum space, it contributes to the matching kernel in the physical region as: ) where F[iπsign(z)] y−x = (y−x) −1 has been used.The real part contribution is proportional to δ(x−y) when convoluted with the jet function, or equivalently, acting as an multiplicative factor to the jet function.In the threshold factorization Eq. (4.6), all the external-state and Dirac-structure dependence has been absorbed into the Sudakov factor.So the jet function has the same form as the quasi-PDF [48], except that the variable is |y − x| instead of |1 − x/y| because the integral measure differs by a factor of 1/|y| in the two cases, and in coordinate space, where l z = ln e 2γ E µ 2 z 2 /4 .In the hybrid scheme, the threshold limit of the momentumspace formalism is modified by [84]: We can verify that the combined contribution reproduces the ERBL logarithm as well as the finite terms of the NLO matching kernel Eq. (4.2) in the threshold limit.
The physical scales in the matching kernel become manifest after threshold factorization.The Sudakov factor incorporates the logarithms of the quark and anti-quark's momenta, and the jet function depends on the soft gluon momentum.Moreover, the quark and antiquark momentum exist in two separate heavy-light Sudakov factors C ± ϕ (p z , µ), thus they can be resummed independently.So in the threshold limit, we are able to resum all the different logarithms in the matching kernel by solving the RG equations.

Resummation of logarithms in the threshold limit
The three parts in the threshold limit are resummed separately.The Sudakov factor follows the RG equation where Γ cusp is the universal cusp anomalous dimension [85] known to four-loop order [86,87], here we use the three-loop result, where γ c is the anomalous dimension of the Sudakov factor known to 2-loop order [48,88,89], The jet function J(∆ = |x − y|P z , µ) evolves as, where γ J is known at up to two-loop order [48], )n f .(4.17) From the RG equation, we can resum the Sudakov factor H(x, P z , µ) numerically by setting an initial scale, such as µ 1 = 2xP z and µ 2 = 2xP z , in C ± (xP z , µ 1 ) and C ∓ (xP z , µ 2 ), then evolve to µ = 2 GeV through numerically solving the differential equation Eq. (4.14).
There is also an analytical solution to the RG equation.The Sudakov factor can be decomposed into the norm and the phase, where the phase is read from the fixed-order expression of and its evolution is just the imaginary part of Eq. (4.14), which can be easily solved, Then the overall phase Â is which turns out to be RG invariant, as the evolution of the imaginary parts in C ± (µ, p z ) and C ∓ (µ, p ′ z ) cancel.So we can write the phase factor in an RG-invariant resummed form in terms of Eq. (4.21), where µ 1 and µ 2 are the physical scales in C ± (µ, p z ) and C ∓ (µ, p ′ z ), respectively.The RGE solution to the norm of the Sudakov factor is [90] where the evolution factors are calculated from the QCD beta function and the anomalous dimensions, The evolution of jet function is more complicated and difficult to implement numerically directly in momentum space.In coordinate space, it is multiplicative, which can be used to derive an analytical solution in momentum space [81], where η = 2a Γ (µ i , µ), and the star function is defined as where ⌊−η⌋ = n for n ≤ −η < n + 1.For η > 0, which happens when µ i > µ, the star function is trivially the function itself.For −1 < η < 0, which is almost always the case for µ i < µ, the star function is the same as the plus function in the range [−∞, ∞].
To implement the threshold resummation, we need first identify the semi-hard scale µ i in Eq. (4.28).It seems natural to set µ i = 2|∆| = 2|y − x|P z , depending on the soft gluon momentum.However, such a choice is dependent on the integrated variable y, and hits the Landau pole at ∆ → 0 for any x value, thus is not numerically implementable.In analogy to the argument in deep inelastic scattering (DIS) [81], we can examine the scale by convoluting the threshold-limit kernel with a DA function and see how the logarithms in the matching kernel are converted to a logarithm depending on the light-cone-DA quark momentum fraction x, after the variable y is integrated out.That is, to examine the following integral when x → 0 and x → 1.If the functional form of ϕ(y) is known, one could perform the integral explicitly, obtaining an analytical function of η, then figuring out the additional logarithm structure generated by the ∂ η operator.It is indeed the case for DIS, when the integral range is limited to y ∈ [x, 1], where the PDF follows a simple power law f (y) ∝ (1 − y) b in the x → 1 limit.The case for DA calculation is more complicated than the PDF case, because the convolution of the kernel with a DA ϕ(y) is always integrating over the full y ∈ [0, 1] range, as shown in Eq. (4.31), and it does not follow the simple power-law form in mid-x region.Therefore, the same argument does not hold for DA.Nevertheless, we can argue that the threshold resummation effect is dominated by the integral region y ∼ x in Eq. (4.31).Firstly, the convolution of the jet function and the pion DA in Eq. (4.31) is mostly enhanced when y → x.Meanwhile, if we perform the threshold resummation in the Mellin moment space, it is the ln 2 N and ln N ln µ terms that are resummed in the Wilson coefficients C N →∞,M (z, µ).The anomalous dimensions for C N M (z, µ) in the N → ∞ limit have the following asymptotic behavior indicating that the diagonal ERBL evolution for ⟨ξ N ⟩ is mixed with the threshold logarithm ln N , thus resulting in the effective scale ∼ N z −1 [46], while the off-diagonal parts only contains ln(zµ), and the physical scale is just z −1 .Therefore, the threshold resummation of quasi-DA matching will only change the diagonal Wilson coefficients C N N , which only affects the higher moments ⟨(2x − 1) N →∞ ⟩ of the light-cone DA, determined by the distribution near the endpoint regions x → 0 and x → 1.Moreover, the contribution from the further endpoint, such as from y → 1 to x → 0, is suppressed compared to the y → 0 region, because the corresponding process requires a exchange of hard gluon.Based on the above argument, we can limit our integration domain to be around the nearest endpoint to check the threshold resummation effect.For example, we can integrate over y ∈ [0, 2x] for x → 0, where ϕ(y) can be modeled by some power-law function ϕ(y) ∼ y a , φ(x) where ln x comes from the commutator [∂ η , x a+η ] and is combined with ln µ 2 /(4P 2 z ) to generate the actual physical scale 2xP z .Similarly, when x → 1, we find So a reasonable choice for the semi-hard scale is 2xP z for x < 1/2 and 2(1−x)P z for x > 1/2, i.e., µ i = 2min(x, x)P z .To obtain the threshold-resummed matching kernel, we perform a two-step convolution: we first remove the fixed-order singular threshold terms (i.e., the Sudakov factor and jet functions in the form of δ(y − x) and 1 y−x ) from the fixed-order matching kernel by convoluting an inverse threshold terms with the original fixed-order matching kernel, then add back the corresponding resummed terms by convoluting with the resummed singular threshold terms, where JH = H ⊗ J or J ⊗ H, both reproducing the same singular threshold terms because the factorization is multiplicative in coordinate space.We average them and consider the difference between the two choices as a systematic error in our final results.Here we show the example of JH = J ⊗ H, where the momentum in H labels the parton momentum in the light-cone DA.To extract the light-cone DA from quasi-DA, we need the inverse form of Eq. (4.35), where for a specific momentum fraction x in ϕ(x), the scales are chosen as, Note that we have resummed the renormalon in our perturbative matching kernel C −1 NLO (µ) with LRR to remove the linear correction.The same linear renormalon from the linearly divergent Wilson line self energy now exist in the jet function.To resum the large logarithms in the LRR-improved matching kernel, we then need to resum the renormalon in the jet function in exactly the same way.To implement it to the resummed form Eq. (4.28), it is more straightforward to modify the coordinate space expression Jz at scale µ 0 = 2z −1 e −γ E , then evolve to the scale µ, resulting in a coefficient (1 (4.38)Then the regularized renormalons will cancel between the fixed-order jet function J NLO (µ h ) and the resummed jet function J −1 TR (µ i , µ h ).To turn off the resummation of either the Sudakov factor or the jet function, we can set the corresponding scale to the same as µ in Eq. (4.36).So we can examine the effect of the two resummations independently, as shown in Fig. 12.Here we truncate the resummed results for x ∈ [0.2, 0.8] to avoid α s (2xP z ) or α s (2xP z ) becoming non-perturbative.As one can see, the resummations of the Sudakov factor and jet function have opposite effects on the final x-dependence of the light-cone DA.The higher-order large logarithms of soft quark (antiquark) momentum tend to enhance the endpoint regions of the DA (or equivalently, suppresses the endpoint contributions to physical observables as a convolution of hard kernel and the light-cone DA), while the higher-order large logarithms of soft gluon momentum tend to suppress the endpoint regions.In total, the DA is enhanced at endpoints after threshold resummation when compared to a fixed-order calculation.
In order to check the convergence of perturbation theory after resummation, we vary the hard scale µ h and semi-hard scale µ i by a factor of c = {1/ √ 2, √ 2}.When the contribution from large logarithmic terms at higher orders of perturbation theory are significant, the resummed results will become sensitive to the choice of initial scales, indicating that the perturbation theory does not work or converge, so the matched DA is no longer reliable.We show the scale variation for the hard scale µ h and semi-hard scale µ i in Fig. 13.It is clear from the plot that the results are not very sensitive to the choice of both scales near x = 0.5, but when approaching the endpoints, the scale dependence becomes very large, indicating that the perturbation theory breaks down.From our data, we have the .Scale variation for the hard scale µ h (magenta) and semi-hard scale µ i (green).The uncertainty is small near x = 0.5, but becomes large when approaching the endpoint regions.
largest pion momentum of P z ≈ 1.85 GeV.The corresponding reliable range of prediction is estimated to be x ∈ [0.25, 0.75].To extend the range of the LaMET prediction, we have to go to higher hadron momentum.Since the range is determined by xP z and xP z , with P z = 3.0 GeV data we could be able to reach x min ≈ 0.15.To show that the reliable range does depend on the hadron momentum, we compare the results with P z ≈ 1.6 GeV in Fig. 14.With a smaller momentum, the scale variation grows faster, thus reduces the reliable range of data, suggesting a similar size of x min P z .We also notice that the results in x ∈ [0.3, 0.7] are consistent for the two momenta, suggesting that the power correction is well-controlled here.
In our final result, we consider the systematic error from the choice of z s ∈ [2a, 3a] in the hybrid renormalization scheme, the choice of λ tail = {8.the hard scale and semi-hard scale µ i,h → cµ i,h .For each of these choices, we obtain ϕ i (x) with a statistical error band.The systematic error band is obtained by requiring the total error band to cover all ϕ i (x) results.Figure 15 shows the relative uncertainties as a function of x.We find that the statistical error at x = 0.5 is δ stat ϕ(0.5) ≈ 4%, when including the scale variation, the error increases to δ stat+scale ϕ(0.5) ≈ 6%, and the total systematic error when further including the coordinate-space extrapolation and the difference choices of z s is δ stat+sys ϕ(0.5) ≈ 8%.When approaching the endpoints, the results are dominated by the scale variations, indicating that the perturbation theory becomes unreliable in this region.Noticing the rapid increase of statistical uncertainty when approaching the endpoint, we truncate the result when the systematic uncertainty increases to 10%.This corresponding to a range of x ∈ [0.25, 0.75], within which LaMET is able to make reliable prediction with controllable systematics.The same procedure is applied to O 3 = γ z γ 5 as well.We show the two results and their comparison in Fig. 16.Both results suggest a flat DA in the mid-x region.Comparing the results of two different operators, we find they are consistent within errors.The flatness of our pion DA result with ϕ(0.5) ∼ 1.1 turns out to be qualitatively similar to the prediction from the lightcone sum rule method [91][92][93].

Comparison with data from HISQ ensembles
We compare our final results on DWF ensembles with a similar calculation based on HISQ ensembles at physical quark masses [19,94] in Fig. 17.The lattice spacing a HISQ = 0.076 fm is similar to that in this work, and the momentum P z = 1.78 GeV is close to the largest one P z = 1.85 GeV in this calculation.The DWF results are slightly flatter when compared to the HISQ results, and the distribution at x = 0.5 is ϕ DWF (0.5) = 1.07 (9), lower than the HISQ result ϕ HISQ (0.5) = 1.19 (8).We also take a ratio to the central value of the DWF DA to examine the relative deviation.The momentum-space DA on DWF ensembles is suggesting a slightly lower distribution than the HISQ ensemble within 2σ deviation.Thus, our calculation suggests that the explicit chiral symmetry breaking effect on the shape of pion DA is potential but not significant.To further distinguish different calculations, it is important to reduce the systematic uncertainty in mid-x region, which could be achieved by a more precise measurement of the long-range correlations to constrain the long-tail extrapolation.Meanwhile, this observation is based on the calculation from one lattice spacing, thus the discretization effects is not well-controlled.A more robust conclusion could be drawn after calculating on multiple lattice spacings and extrapolating the results to the continuum.

Conclusion
In this article, we present the first lattice QCD calculation of the x-dependence of pion DA with the chirally symmetric DWF fermion action at physical pion mass.For the first time, we reformulate the threshold resummation for the quasi-DA matching kernel, and resum both the ERBL and threshold logarithms with this technique in the perturbative matching to extract the light-cone DA.We analyze two different Dirac structures γ t γ 5 and γ z γ 5 , both of which approach the light-cone limit γ + γ 5 when boosted to infinite momentum P z → ∞.With the largest pion momentum P z ≈ 1.85 GeV, we demonstrate that the region x ∈ [0.25, 0.75] is not sensitive to scale variations, where the systematics is under control.Beyond this region, the scale variation becomes much larger, indicating that the perturbative matching is no longer reliable.We also show that increasing the pion momentum could extend the reliable range of x.Our results are consistent between the two different Dirac structures.Applying the same analysis method to the quasi-DA from a HISQ ensemble at physical point, we notice that the pion DA from the DWF ensemble is flatter, with a smaller amplitude around x = 0.5.The discrepancy is within 2σ, thus is not statistically significant enough to be conclusive.It should be pointed out that our study is limited to one lattice spacing, which cannot estimate the discretization effects until a continuum extrapolation is performed.Extending the future work to finer lattice spacings to examine the results in the continuum limit will give us a more convincing and comprehensive understanding of the chiral symmetry breaking effects on pion structure.
A z s -dependence in the hybrid renormalization scheme The hybrid scheme introduces a piecewise function Eq. (3.11), that guarantees the continuity but not smoothness of the renormalized matrix elements.Also, it introduces z s -dependence to the calculation.This is not a problem in practice, because the perturbative matching kernel in Eq. (4.2) also depends on z s , and could cancel the z s -dependence in quasi-DA, as shown in the left plot of Fig. 18.
Meanwhile, we can also easily connect the renormalized matrix elements in hybrid schemes with different z s values through a perturbative correction, assuming z ′ s > z s > 0: where C 0 (z, µ) is the Wilson coefficient evaluated at z.Then, after this perturbative correction, we find a nice agreement for z s = 2a and z s = 3a in coordinate space, as shown in Fig. 18.When z s is small, the C 0 (z, µ) depends logarithmically on z s .Thus we would expect that the "non-smoothness" in the coordinate-space renormalized matrix elements comes from ∂C 0 (z, µ)/∂z s , which would be suppressed by z.When we choose larger z s , the curve will look more smooth in Fig. (10).But as we have just shown above, the choice of z s doesn't really affect our calculation.

9 Figure 1 .
Figure 1.Effective mass from the C ππ (t, 0) correlators.The straight lines represent the energy calculated from dispersion relation E(P z ) = P 2 z + m 2 .

Figure 2 .
Figure 2. First excited state energy from two-state fit as function of t min for different pion momenta.

Figure 3 .
Figure 3. Examples of fitted correlators at z = 0 (left) and z = 3a (right) with different t min .

Figure 4 .
Figure 4. Dispersion relation for ground state energy E 0 and the first excited state energy E 1 , fitted with priors.The blue and orange bands are the prior for E 0 and E 1 .

9 Figure 5 .
Figure 5. Bare f π obtained from two-state fits for pion with different momenta as function of t min .

Figure 6 .
Figure 6.Matrix elements from two-state fit of C πO0 (top) and C πO3 (bottom).

Figure 7 .
Figure 7.The fitted ratio (left) and the second moments (right) from the fits at different z values.

Figure 9 .
Figure 9. Renormalized matrix element (left) and the RG-invariant ratio (right) compared with OPE reconstruction in dashed lines at short distance for O 0 (top) and O 3 (bottom).Different lines correspond to different momenta.When momentum increases, the renormalized matrix elements or the ratio goes down.

Figure 10 .
Figure 10.Real part of renormalized matrix elements for O 0 (left) and O 3 (right).At large zP z , the data from different momenta follow a general curve, except for scaling violation in z 2 or 1/P 2 z

Figure 11 .
Figure 11.long-tail extrapolation and the corresponding x-dependent quasi-DA.The different long-tail models introduce up to 10% systematic uncertainties in the x dependence.

Figure 12 .
Figure 12.Matching from quasi-DA to light-cone DA with fixed order inverse kernel (black), resummed kernel (blue), and partially resummed kernels for the Sudakov factor (magenta) and the jet function (green).
Figure 13.Scale variation for the hard scale µ h (magenta) and semi-hard scale µ i (green).The uncertainty is small near x = 0.5, but becomes large when approaching the endpoint regions.

5 ,Figure 14 .
Figure 14.Left: hard scale variation for smaller pion momentum P z = 1.6 GeV.Right: a comparison between results from 1.8 GeV 1.6 GeV.

Figure 15 .
Figure 15.Relative uncertainties as a function of x.

Figure 16 .
Figure 16.Final light-cone DA for O 0 (left) and O 3 (middle), including scale variation and systematic errors.And a comparison between the two operators are shown on the right.Outside the range x ∈ [0.25, 0.75], we cannot make a reliable prediction with perturbative matching.

Figure 17 .
Figure 17.Comparison between the pion DAs calculated on DWF and HISQ ensembles (left) and their ratios to the central value of the DWF DA (right).

Figure 18 .
Figure18.The z s dependence of hybrid scheme in momentum space (left) and coordinate space (right).Although the quasi-DA are different with different z s choices.The momentum-space perturbative matching or a coordinate-space perturbative correction could cancel such z s dependence.