Pion Valence-Quark Generalized Parton Distribution at Physical Pion Mass

We present the first lattice-QCD $x$-dependent pion valence-quark generalized parton distribution (GPD) calculated directly at physical pion mass using the Large-Momentum Effective Theory (LaMET) with next-to-next-to-leading order perturbative matching correction. We use clover fermions for the valence action on $2+1+1$ flavors of highly improved staggered quarks (HISQ), generated by MILC Collaboration, with lattice spacing $a \approx 0.09$fm and box size $L \approx 5.5$fm; the pion two-point measurements number up to $O(10^6)$ with boost momentum 1.73GeV. The pion valence distribution is renormalized in hybrid scheme with Wilson-line mass subtraction at large distances in coordinate space, followed by a procedure to match it to the $\overline{\text{MS}}$ scheme. We focus on the zero-skewness limit, where the GPD has a probability-density interpretation in the longitudinal Bjorken $x$ and the transverse impact-parameter distributions. We take the integral of our GPD functions to generate leading moment so that we can make comparisons with past lattice-QCD and experimental determinations of the pion form factors and found consistent agreement among them. We predict the higher GPD moments and reveal $x$-dependent tomography of the pion for the first time using lattice QCD.


I. INTRODUCTION
Pions, as the Nambu-Goldstone bosons created by the chiral symmetry breaking of quantum chromodynamics (QCD), play a crucial role in our understanding of the origin of mass for matter [1][2][3][4].Decades of effort have been devoted to pursuing understanding of how quarks and gluons give rise to pions.Among the various properties of the pions, generalized parton distributions (GPDs), a concept introduced more than two decades ago [5,6], have been among the most deeply studied both theoretically and experimentally.(We refer interested readers to the review papers [7][8][9][10].)GPDs retain information on the form factors' dependence on the transfer momentum Q 2 as well as the one-dimensional structure in terms of the Bjorken-x parameter in the parton distribution functions (PDFs).Multiple exclusive processes, such as deeply virtual Compton scattering (DVCS) and deeply virtual meson production (DVMP), provide experimental access to the GPDs, allowing us to map out the tomography of hadrons [11].DVCS has been extensively studied during the last two decades; the process is relatively simple in its interpretation, since the only nonperturbative object entering its amplitudes, which are referred to as Compton form factors (CFFs), are the GPDs.They are the main source of information on GPDs.We refer readers to Ref. [9] for phenomenological context.Efforts have been made to extract information experimentally about the internal structure of pion, ranging from its electromagnetic form factor through pion-electron scattering [12] to its PDFs through the Drell-Yan (DY) process [13][14][15].As a result, there has been much progress made toward exploring the global analysis of pion PDFs in recent years [16][17][18], but not much effort has been made toward pion GPDs.There are many ongoing and planned experimental efforts to further our knowledge of pion structure; for example the JLab 12-GeV program will provide precision data relating to the pion and kaon form factors up to Q 2 ≈ 10 GeV 2 and 5 GeV 2 , respectively.They will also measure the pion and kaon structure functions at x > 0.5 through the Sullivan process [19].AMBER Collaboration at CERN can play a crucial role due to their unique capability delivering pion and kaon DY measurements in the centre-of-mass (CoM) energy range 10-20 GeV.A future facility, the Electron-Ion Collider (EIC) at Brookhaven National Laboratory, will provide information on pion and kaon structure over a large and tunable CoM energy range 20-140 GeV via Sullivan-process measurements [1,2,4,10,20,21].An electron-ion collider being discussed in China (EicC) [22,23] would have a similar CoM energy range to AMBER and neatly fill a gap between JLab 12 and the EIC.
Less progress has been made toward x-dependent GPD studies by contrast.The first lattice study of GPDs was made for pions [45] in 2019 with largest boost momentum 1.7 GeV using 310-MeV pion mass.By 2020, ETMC [57] reported results at a single Q 2 for both unpolarized and polarized GPDs with a 260-MeV pion with skewness of 0 and 0.3 with a single transfer momentum of 1.67 GeV.In the same and following year, MSULat [58,61] reported on the unpolarized and polarized zero-skewness GPDs with 135-MeV pion at multiple Q 2 .One can actually take the moments of MSULat's x-dependent GPDs and compare with lattice results using a traditional moment calculation using a local operator via the OPE method, such as electromagnetic and axial form factors (with n = 1), and generalized form factors (n = 2).With multiple Q 2 values in the zero-skewness limit, the Fourier transform of the non-spin-flip GPD H(x, ξ = 0, −Q 2 ) yields the impact-parameter-dependent distribution, which cannot be directly accessed via experimental measurements.Progress so far has been done using Breit frame, where the initial and final momenta of the hadron differ by half the transfer momentum.Recent work on asymmetric-momentum setups for GPDs has been demonstrated by the ETM and BNL/ANL groups [62] and can help reduce the computational cost of the lattice calculation.Pion structure, such as form factors, can be more sensitive to the pion mass used in the lattice calculation; therefore, it is important to study pion structure directly at the physical pion mass.
In this work, we report the first study at physical pion mass of the pion unpolarized valence-quark GPD, using LaMET method with lattice spacing 0.09 fm in the Breit frame.The remainder of this paper is organized as follows: In Sec.II, we discuss the lattice setup and the procedure for how the lattice two-and three-point correlators are calculated and analyzed to extract the ground-state matrix elements.These matrix elements are then renormalized using hybrid renormalization scheme and the valence-quark GPDs are extracted by fitting the matrix elements with NNLO kernels in Sec.III.We compare our results in certain limits that have prior LQCD calculations for comparison: our valence-quark GPD at Q 2 result is consistent with prior lattice physical-pion-mass calculate the quasi-GPD distribution and matched them to the lightcone, and the first moment of our GPD is in good agreement with the prior lattice pion form factors.We predict the transfer-momentum dependence of the pion GPDs, higher moments of the generalized form factors that have not been calculated on the lattice, and the tomography of the pion.Conclusions and future outlook can be found in Sec.IV.

II. CALCULATION SETUP
The unpolarized valence-quark GPD of the pion on the lightcone is defined as where Bjorken-x is the momentum fraction x ∈ [−1, 1], µ is the renormalization scale in the MS scheme, momentum P µ = (P 0 , 0, 0, P z ), where the pion momentum of the initial and final states are P ∓ ∆/2 with ∆ the momentum transfer, and the variables t = ∆ 2 and ξ = − ∆ + 2P + .q and q represent the antiquark and quark fields, and the gauge link ) ensures gauge invariance of the quark bilinear operator.
In the forward limit (∆ µ → 0), it reduces to the PDF.
On the lattice, one can compute ground-state matrix elements of a pion with a finite-P z boost in the Breit frame, where U z is a discrete gauge link in the z direction, ⃗ P = {0, 0, P z } is the 3-momentum of the pion, Γ = γ t and ⃗ ∆ is the momentum transfer between initial and final pion.The matrix elements are renormalized and then contribute to the quark quasi-GPD via Hπ q (x, ξ, t, P z , μ) = using the renormalized matrix element hR .The lightcone GPD in the MS scheme at scale µ is then convolved with a perturbative hard matching kernel, up to power corrections that are suppressed by the pion momentum [116] Hπ two-sim FIG. 1.
An example of the bare ground-state matrix-element determination using two-state fits to multiple source-sink separations, along with the ratio plot and reconstructed fit using The colored bands are the reconstructed ratios from the fits to source-sink separation tsep ≈ [0.72, 0.9] fm, respectively, and the gray band shows the ground-state matrix elements from the "two-simRR" fits.
In the zero skewness limit ξ = 0, the matching kernel C is the same as the matching kernel for the PDF, as documented in Refs.[37,117].
This calculation is carried out using N f = 2 + 1 + 1 highly improved staggered quarks (HISQ) [118], generated by the MILC Collaboration [119], with lattice spacing a ≈ 0.09 fm at physical pion mass.Five steps of hypercubic (HYP) smearing [120] are applied to the gauge links to reduce short-distance noise.We use Wilson-clover fermions in the valence sector and tune the valence-quark masses to reproduce the lightest light and strange sea pseudoscalar meson masses.A similar setup is used by PNDME collaboration [121][122][123][124][125][126][127][128][130][131][132] with local operators, such as isovector and flavor-diagonal charges, form factors and moments; the results from this mixed-action setup are consistent with the same physical quantities calculated using different fermion actions [60,[133][134][135][136][137].BNL/ANL Collaboration has also been using a mixed-action setup for their pion valence-quark PDFs calculations [110,114].We carefully monitor the measurements from each configuration to rule out any possibility of exceptional configurations; we have not observed any here nor in the previous works mentioned earlier.
To calculate the GPD matrix elements at nonzero momentum transfer, we first calculate the matrix element ⟨χ π ( ⃗ P f )|O µ |χ π ( ⃗ P i )⟩, where χ π is the pion interpolating field, O µ = ψγ µ W (z)ψ is the LaMET Wilson-line displacement operator with ψ being the quark field, and ⃗ P {i,f } the initial and final nucleon momenta.We use 1960 configurations to help with statistical noise at boost momentum L {0, 0, 8}a −1 .We vary the spatial momentum transfer ⃗ q = ⃗ P f − ⃗ P i = 2π L {n x , n y , 0}a −1 with integer n x,y and n 2 x + n 2 y ∈ {0, 4, 8, 16, 20} with four-momentum transfer squared t = Q 2 = −q µ q µ ∈ {0, 0.19, 0.39, 0.77, 0.97} GeV 2 using periodic boundary conditions.We use source-sink separations of {6,7,8,9,10} lattice units with high-statistics measurements of {752,640, 1,003,520, 1,003,520, 1,505,280, 1,505,280}, respectively.We extract pion ground-state matrix elements using two-state simultaneous fits about which more details can be found in our previous work [40,58,61,126].Figure 1 shows an example of the ground-state matrix elements (shown as the gray band) from the simultaneous two-state fitted results using source-sink separation of t sep ∈ [6,10] in lattice units, at L .The signal in the three-point correlators decreases as the source-sink separation increases; also as expected, the ground state for each t sep approaches the simultaneous-fit ground-state matrix element, shown as the gray band in the plot.The spectral decomposition predicts that the data for all three quantities is symmetric about t = t sep /2 in Breit frame; at large t sep such symmetry breaks slightly due to the statistical fluctuations.The groundstate matrix elements are consistent when one removes smaller t sep from the analysis.

III. RESULTS AND DISCUSSION
We use the hybrid renormalization scheme [138] with NNLO perturbative matching correction [51,115].The exact procedure has been applied and detailed in prior valence-quark PDF studies by BNL/ANL group [110,114].This renormalization scheme is a combination of the ratio scheme at short distances up to Wilson-line length z s and an explicit subtraction of the self-energy divergence of the Wilson line at large distances.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.For our calculations, we choose z s = 0.27 fm. .The corresponding bands are the reproducing matrix elements from the fit to the matrix elements using Eq. 4. (Right) Comparison of two methods, fit-matching "fit" and direct matching with extrapolation "extrap", to extract the lightcone distributions for valence-quark pion GPD at Q 2 = 0 and 0.97 GeV 2 .The two methods are consistent over most x values within statistical errors.There is some tension near x = 0.25 and x > 0.9 for H π (x, Q 2 = 0) GPD, but they are still consistent within two sigma.
The hybrid normalized matrix elements as functions of zP z can be found on the left-hand side of Fig. 2. We also normalize all the matrix elements by the bare ground-state matrix element at z = 0 and Q 2 = 0 (vector charge normalized to 1) to improve the signal.The real matrix elements decrease quickly to zero as a function z due to the large boost momentum used in this calculation.The signal-to-noise ratio also increases quickly as one increases the Wilson-line displacement.The central values of matrix elements also drop quickly as Q 2 increases initially but slow down at the larger Q 2 ; this is expected, since similar behavior observed in the Q 2 -dependence of the pion form factors.
The valence-quark GPD H π (x, Q 2 ) is then extracted using Eq. 4 with NNLO matching kernel C from Refs.[51,115].We compare two methods for performing the matching.First, we extrapolate the matrix elements to large λ = zP z using a form with the exponential and power decay, d 1 e −d2|z| /|λ| d3 , as used in prior pion-PDF studies [110,114].The fit parameters d 1 , d 2 and d 3 vary with Q 2 and are fit to the data subset at large zP z > 3.9.We then Fourier transform the matrix elements into the quasi-distribution and apply matching according to Eqs. 3 and 4. Alternatively, we can adopt a phenomenology-inspired functional form to describe the x dependence of the lightcone GPD: where the parameters m, n and c are Q 2 -dependent, and the beta function B(m + 1, n + 1) = 1 0 dx x m (1 − x) n is used to normalize to unity at Q 2 = 0. Note that this form is also commonly used by the PDF global-fit community, for example, in the pion PDF fit by JAM [16,139].This phenomenological form has also been widely used for xdependent valence-quark pion PDFs, such as in Refs.[49,86,87,95,110,140].We perform the phenomenological fit to our physical pion-mass ensemble using the matrix elements with Wilson-line displacement z ∈ [0.09, 0.81] fm for each transfer momentum.We choose the MS renormalization scale to be µ = 2 GeV.We find the goodness of fit to be around 1 with either c as a free parameter or c = 0.As in the previous lattice valence-quark pion PDF study, the final distributions are consistent with either choice, but the former case results in the distribution having larger error bands due to the additional free parameter.For the remainder of this work, we focus on the c = 0 results, since this is the first lattice work on the pion GPD at physical pion mass.Further study of the systematics due to fit-form choices, as well as other lattice artifacts (lattice discretization, large momentum dependent, etc.), should be explored when more computational resources are available.The reconstructions of the matrix elements using the fit parameters are shown as the bands on the left-hand side of Fig. 2; they pass through the input data nicely.Varying z max results in small changes to the goodness-of-fit and x dependence of the distributions.The two methods of extracting the GPD are consistent for all Q 2 values studied here.We show two selected values of Q 2 results on the right-hand side of Fig. 2, labeled as "extrap" and "fit" for the former and latter methods, respectively.Both results are consistent in the intermediate-x region within statistical errors.Some small disagreement appears for x around 0.25 and at large x of H π (x, Q 2 = 0) but still within two standard deviations; this may change as the statistics increase in future updates.For the remainder of the manuscript, we focus on the results from the fitting method.Figure 3 shows our final result for the pion valence-quark GPD.On the left-hand side of Fig. 3, the GPD at zero transfer momentum as a function of Bjorken-x (that is, the PDF) is shown, as well as the distribution at various transfer momenta in the zero-skewness limit.As the transfer momentum increases, the distribution decreases, as expected.On the right-hand side of Fig. 3, we show the GPD at selected x values as functions of transfer momentum Q 2 (in units of GeV 2 ).The GPD is fitted with a z-expansion up to 3 parameters to interpolate the Q 2 dependence.The Q 2 dependence of the GPD can help us investigate the tomography of the pion.Unfortunately, there are few published lattice pion GPDs at physical pion mass to make comparisons.In first lattice GPD work [45], the pion valence-quark GPD at zero skewness was calculated using clover valence fermions on an ensemble of gauge configurations with lattice spacing a ≈ 0.12 fm, box size L ≈ 3 fm and pion mass m π ≈ 310 MeV.Our results are consistent with the heavier quark mass but with the systematics removed or improved for being at the physical pion mass and smaller lattice spacing.We can compare our results at zero transfer momentum with the PDF results of BNL/ANL group, also done at physical pion mass but at a finer lattice spacing of 0.076 fm [110,114]; they are consistent within statistical errors.
and pion masses ranging from 139 to 340 MeV.Their statistics range from 9,600 to 485,376 measurements [154].BNL used clover on N f = 2 + 1 HISQ lattice at a single lattice spacing a = 0.076 fm at physical pion mass with 1750 and 35,000 exact and sloppy inversions, respectively [155]; the work also includes two additional ensembles at smaller lattice spacing with 300-MeV pion to quantify the systematics.They found the pion radius to be consistent at heavy pion mass between the 0.04 and 0.06 fm lattice-spacing ensembles but sensitive to the pion mass.We select their results at physical pion mass on the smallest lattice-spacing to compare to our integral over the GPD.We can see that there is a very nice agreement among all the lattice results, even though the lattice spacing varies among the three calculations.We also compare the pion form factor with those extracted from experiments [157][158][159][160][161] and find good agreement; the lattice data are more precise than the available experimental data in certain transfer-momentum regions.
There have been only a couple of dynamical calculations of A π 2,0 using heavy pion masses more than a decade ago, reported in conference proceedings.In 2005, UKQCD/QCDSF [162] used two flavors of dynamical fermion with a very heavy pion mass of 1090 MeV and found A π 20 (t = 0) = 0.315 (8).In 2013, RQCD also used two-flavor clover dynamical fermions with pion masses ranging from 150 to 491 MeV [163]; unfortunately, they found A π 20 (t) to have a wide range distribution and the data at the near-physical pion mass has uncertainty too large for a direct comparison.However, overall, our moment-integral from GPD functions lies within the range of RQCD moment results obtained using OPE method.Unfortunately, there have been no n ≥ 3 GFF moment calculations from OPE methods.Our A π {3,4}0 results, shown in Fig. 4, are the only predictions from LQCD.
Taking the lattice calculations of the pion's valence quark GPD, H π (x, ξ = 0, Q 2 ), we can then Fourier transform of this GPD to learn about the impact-parameter-dependent distribution, q(x, b) [7,165]: where b is the light-front distance from the center of transverse momentum (CoTM).Figure 5 shows slices of the distribution with selected x ∈ [0.1, 0.9], as well as the two-dimensional distribution at x = 0.45.The impactparameter-dependent distribution describes the probability density for a parton with momentum fraction x to be found in the transverse plane at distance b from the CoTM.It provides a snapshot of the pion in the transverse plane and indicates what might be expected from pion tomography.

IV. SUMMARY AND OUTLOOK
We presented a state-of-the-art high-statistics lattice calculation of the valence-quark GPD of the pion using the LaMET approach with a next-to-next-to-leading order perturbative matching formula.The calculation was done at lattice spacing 0.09 fm with physical pion masses and boosted pion momentum around 1.7 GeV with four additional nonzero transferred momenta in the Breit frame.We performed multistate analyses with multiple source-sink separations to remove excited-state contamination.We obtained the x-dependent valence-quark GPD by fitting the hybrid-scheme-renormalized matrix elements to a phenomenology-inspired functional form.We compared our results with prior lattice studies in special limits of the pion GPD, such as the pion PDFs (Q 2 = 0 limit) and pion form factors (forward limit).This is the first pion valence-quark GPD directly calculated at the physical quark masses.We also made predictions for higher moments of the generalized form factors that have not yet been calculated and show x-dependent tomography of the pion for the first time using lattice QCD.Our result provides the most reliable lattice prediction of the valence-quark pion GPD to date and will guide upcoming measurements at JLab and the EIC.Future work will extend the calculation to additional boost momenta and finer lattice spacings to further constrain the systematics in the lattice calculations.

Q 2 (GeV 2 ) 5 Q 2 (FIG. 3 .
FIG. 3.(left) Pion valence-quark GPD as a function of Bjorken-x with five values of transfer momenta studied in this work.(right) Pion valence-quark GPD as a function of transfer momenta at selected Bjorken-x indicated in the bands.The z-expansion is used to interpolate between the five transfer momenta.

45 FIG. 5 .
FIG. 5. (left) The valence-quark impact-parameter-dependent distribution of pion as a function of b for selected x and (right) an example of the two-dimensional distribution at x = 0.45.All distances are in units of fm.
FIG. 2. (Left)The hybrid renormalized Wilson-line length displacement z-dependence of matrix elements as a function of dimensionless parameter zPz at various momentum transfer Q 2 ∈ {0, 0.97} GeV 2