Unpolarized Transverse-Momentum-Dependent Parton Distributions of the Nucleon from Lattice QCD

a first lattice QCD calculation of the unpolarized nucleon’s isovector transverse-momentum-dependent parton distribution functions (TMDPDFs), which are essential to predict observables of multi-scale, semi-inclusive processes in the standard model. We use a N f = 2 + 1 + 1 MILC ensemble with valence clover fermions on a highly improved staggered quark (HISQ) sea to compute the quark momentum distributions in a large-momentum nucleon on the lattice. The state-of-the-art techniques in renormalization and extrapolation in the correlation distance on the lattice are adopted. The perturbative kernel up to next-to-next-to-leading order is taken into account, and the dependence on the pion mass and the hadron momentum is explored. Our results are qualitatively comparable with phenomenological TMDPDFs, which provide an opportunity to predict high energy scatterings from first principles.


INTRODUCTION
Since the nucleon is at the core of atoms and accounts for nearly all of the mass of the visible universe, exploring its internal structure has been a key task for more than a century in both particle and nuclear physics.In high-energy scattering, the quark and gluon transverse momentum and polarization degrees of freedom in the nucleon are best described by transverse-momentum parton distribution functions (TMDPDFs).Thus, mapping out the nucleon's TMDPDFs is a crucial step in understanding the interactions between quarks and gluons, and possibly the phenomenon of confinement [1,2].Moreover, predicting the observables in multi-scale, noninclusive high energy processes such as semi-inclusive deep-inelastic scattering and Drell-Yan scattering at the large hadron collider (LHC) or electron ion collider (EIC) heavily relies on the knowledge of TMDPDFs [3,4].
Whereas high energy experiments have accumulated a wealth of relevant data, our knowledge of TMDPDFs is far from being complete.Their rapidity evolution, i.e. the Collins-Soper kernel [1], has been perturabtively calculated up to four loops [5,6], but TMDPDFs at low energies are nonperturbative in nature.Based on thousands of data points from the low-p T semi-inclusive DIS and Drell-Yan scattering processes and perturbative-QCD factorization, a number of phenomenological analyses have been made to obtain state-of-art TMDPDFs [7][8][9][10][11].While similar datasets were employed in these analyses, the outcomes exhibit notable discrepancies.This suggests the presence of significant uncertainties in the global extraction of TMDPDFs, underscoring the need for additional constraints to achieve a more refined determination.

arXiv:2211.02340v3 [hep-lat] 19 May 2024
First-principles calculations of TMDPDFs require nonperturbative methods such as lattice QCD.A handful of available investigations using lattice QCD are limited to the ratios of moments of TMDPDFs [12][13][14][15].The development of large momentum effective theory (LaMET) allows the extraction of light-cone quantities through the simulation of equal-time quasi distributions [16,17].A directly caluclation is TMDPDFs is non-trivial due to the presence of the soft function [18], which involves two opposite light-like directions.Implementing this on an Euclidean lattice is a crucial difficulty.Recent progress demonstrates that the rapidity-independent (intrinsic) soft function can be calculated from a large-momentumtransfer form factor of a light meson [19], while the rapidity evolution kernel in the soft function can be accessed via the quasi TMDPDFs/beam functions [18,[20][21][22] or quasi transverse-momentum-dependent wave functions [19,23].Subsequent lattice efforts have been devoted to exploring the Collins-Soper kernel and intrinsic soft function.The agreement between lattice results and phenomenological analyses is encouraging [21,[24][25][26][27].
Following these developments, this work presents a first calculation of TMDPDFs from first principles.We simulate the TMD momentum distributions in a large momentum nucleon or quasi TMDPDFs on the lattice and perform a systematic study of renormalization properties by considering the subtractions from a combination of Wilson loop and short distance hadron matrix element [28].In the matching from quasi TMDPDFs, we include one-loop and two-loop perturbative contributions and employ the renormalization group equation to resum the logarithms.After analyzing the valence pion mass and momentum dependence, our final results for TMDPDFs are found to have a similar behavior as phenomenological fits.
The remainder of this paper is structured as follows.Sec.II presents the theoretical framework, followed by the presentation of Lattice simulations in Sec.III.Sec.IV details the final results for TMDPDFs, while a concise summary and future prospects are outlined in Sec.V.

A. Constructing the equal-time quasi TMDPDFs
Describing the momentum distributions of a parton inside a hadron, TMDPDFs f (x, b ⊥ , µ, ζ) are functions of the longitudinal momentum fraction x, the Fourier conjugate b ⊥ of the parton transverse momentum q ⊥ , as well as the renormalization scale µ and the rapidity scale ζ.In this work we will consider the flavor non-singlet/isovector unpolarized quark TMDPDFs, which do not mix with gluons.
In LaMET, the correlations with modes traveling along the light-cone can be extracted from distributions in a fast-moving nucleon through large-momentum expan-sion.On the lattice, the equal-time quasi TMDPDFs are constructed as where a denotes the lattice spacing.Γ = γ t or γ z is the Dirac matrix that can be projected onto γ + in the large momentum limit.The h0 Γ (z, b ⊥ , P z , a, L) is built with a gauge-invariant nonlocal quark bilinear operator as In the above, |P z ⟩ denotes the unpolarized nucleon state and L denotes the "infinity" that the gauge link can reach.The staple-shaped Wilson link is chosen as in which U z is the path-ordered Euclidean gauge link along the z-direction, and U ⊥ is along the transverse direction at the "infinity" position (z + L) on the finite lattice.Their explicit forms are: An illustration of the quasi TMDPDFs is given in Fig. 1, in which The staple-shaped Wilson link is depicted as double lines.
Illustration of the quasi TMDPDFs.

B. Renormalization of the TMD matrix elements
Quantities in Eq. ( 2) and (3) with the superscript "0" are bare quantities on a finite lattice.They contain linear divergence, pinch-pole singularity and logarithmic divergence.Both the linear divergence and pinch-pole singularity can be renormalized by the square root of Wilson loop Z E (2L + z, b ⊥ , a) [29][30][31][32][33]; and the logarithmic divergence can be renormalized by a factor Z O (1/a, µ), which is extracted from matching between the lattice and perturbative calculation of zero-momentum matrix elements in the perturbative region [28,[34][35][36].
The Wilson loop Z E (r = 2L + z, b ⊥ , a) is defined as the vacuum expectation of a rectangular shaped spacelike gauge links with size r × b ⊥ , which can be written as The Wilson loop is introduced to eliminate the linear divergence of the form e −δ mr , which comes from the selfenergy corrections to the gauge link [29,35], as well as the pinch-pole singularity, which is related to the heavy quark effective potential term e −V (b ⊥ )L describing the interactions between the two Wilson lines along the z direction in the staple link [22].
The logarithmic divergence factor Z O can be extracted from the zero-momentum bare matrix elements h0 Γ (z, b ⊥ , 0, a, L).In order to keep the renormalized matrix elements consistent with perturbation theory, Z O should be determined from the condition: in a specific window where z ≪ Λ −1 QCD so that perturbation theory is valid.The perturbative result for the Γ(= γ t or γ z ) zero-momentum matrix element up to oneloop order in the MS scheme reads Here the perturbation results have been evolved from the intrinsic physical scale µ 0 = 2e −γ E / z 2 + b 2 ⊥ to the MS scale µ using the renormalization group equation [37].
C. Matching at next-to-next-to-leading order and renormalization group resummation It has been shown that quasi TMDPDFs have the same collinear degrees of freedoms as TMDPDFs [22].Their differences from soft modes can be attributed to the intrinsic soft function and different rapidity scales.Also contributions from highly off-shell modes are local [20].Thus TMDPDFs f (x, b ⊥ , µ, ζ) are connected to quasi TMDPDFs fΓ (x, b ⊥ , ζ z , µ) via a multiplicative factorization [22,38]: where S I denotes the intrinsic soft function that has been calculated on the lattice in Refs.[25,26,39,40]  The matching kernel H, as a function of 2 /µ 2 , has been perturbatively determined up to next-to-leading order (NLO) [20,22,41,42] H (1) as well as next-to-next-to-leading order (NNLO), calculated recently [43,44] where ζ z = (2xP z ) 2 and c 2 = 0.0725C 2 F − 0.0840C F C A + 0.1453C F n f /2.The perturbative TMDPDF is calculated in the MS scheme with a fixed renormalization scale µ, while the quasi TMDPDFs are associated with the Collins-Soper scale √ ζ z , which is the intrinsic physical scale of perturbative matching.In order to expose the intrinsic physical scale [37], we resum the large logarithms ∼ ln n ζ z /µ 2 in the small x region through the renormalization group (RG) equation for H: where  3 + 160 27 [22,42].The cusp anomalous dimension Γ cusp is known up to the four-loop level for the quark case [43,45,46].
In practice, we compare the matching kernel at fixed order and with employing the RG resummation starting from the Collins-Soper scale µ 0 = 2xP z to µ = 2GeV.After the resummation, the intrinsic scale 2xP z appears in the running coupling α s (2xP z ).Fig. 2 shows the comparison of NLO and NNLO matching kernels with and without RG evolution.One can see that the RG evolution changes the perturbative behavior at small-x, and makes the predictions of the TMDPDFs in this region less reliable.

III. LATTICE SIMULATIONS
We use the valence tadpole improved clover fermion on the hypercubic (HYP) smeared [47] 2 + 1 + 1 flavors MILC configurations with highly improved staggered quark (HISQ) sea and 1-loop Symanzik improved gauge action [48].We analyze a single ensemble with lattice spacing a = 0.12 fm and volume n 3 s × n t = 48 3 × 64 using physical sea quark masses, and two choices of light valence quark mass corresponding to m val π = {220, 310} MeV.The HYP smearing is also used for nonlocal correlation functions to improve the statistical signal.In order to explore the momentum dependence, we employ three different nucleon momenta as P z = 2π/(n s a) × {8, 10, 12} = {1.72,2.15, 2.58} GeV.
We adopt momentum-smearing point source [49] at several time slices, and average correlation functions for both the forward and backward directions in z and transverse space of the gauge link.In total, there are 1000 (configurations) ×16 (source time slices) ×4 (forward/backward directions of the z and transverse axes) measurements for the m val π = 220 MeV case and 1000 × 4 × 4 measurements for the 310 MeV case.

A. Dispersion relations
The two-point correlation functions of the nucleon are defined as where T = (1 + γ t ) /2 denotes the unpolarized projector, and χ = ϵ abc u a u T b Cγ 5 d c is the nucleon interpolation field.We have computed the two-point functions using nucleon momenta up to 2.58 GeV to examine the dispersion relation.
For the three momenta P z = 1.72, 2.15, 2.52 GeV, the statistics of the two-point correlation functions are 1000 (configurations) ×16 (source time slices) for m π = 220 MeV and 1000 × 4 for m π = 310 MeV; while for the cases with a momentum smaller than 1.72 GeV, there are 1000 × 4 measurements for m π = 220 MeV and 1000 × 1 measurements for m π = 310 MeV.Throughout we use the parametrization C 2 (t) = c 0 e −E0t 1 + c 1 e −∆Et , and we perform two-state fits to extract the ground-state energies, as shown in Fig. 3.
Based on the ground-state energies with different P z , one can obtain the dispersion relations of the nucleons with the pion mass m π = 220 MeV and 310 MeV.We adopt the following parametrization where the quadratic term of lattice spacing a is introduced to parameterize discretization errors.The fit results are shown in Fig. 4. For the m π = 220 MeV case, it is found that b 1 = 1.014(95) and b 2 = −0.014(17), while for the 310 MeV case, the fit gives b 1 = 1.066(80) and b 2 = −0.015(14).From these results, we can see that the dispersion relation is consistent with E (P z ) = m 2 + (P z ) 2 within uncertainties.

B. Bare quasi TMDPDFs from correlated joint fits
To extract the quasi TMDPDFs, one constructs the three-point functions We adopt the sequential source method [50] to reduce the number of propagators in three-point functions, t s denotes the time position of the sequential source.The operator ÕΓ TMD (⃗ x, t) is short for the TMD nonlocal quark bilinear operator Õ0 The three-momentum is chosen as ⃗ P = (0, 0, P z ).
After inserting the single particle intermediate states, one can parametrize the ratio of three-and two-point functions as in which h0 Γ ≡ h0 (z, b ⊥ , P z , Γ), ∆E is the mass gap between the ground-state and excited state, and c 1,2,3 are parameters for the excited-state contaminations.Combining the parametrization form of 2pt functions:

one can extract the values of h0
Γ at fixed (z, b ⊥ , P z ) through a correlated joint fit.
In our calculation, we use bootstrap resampling to establish correlations among all datasets, and the correlations are maintained consistently throughout the entire analysis.For the three-point functions, we have excluded the contact points (t = 0 and t = t seq ) for the five values of the source-sink separation time range t seq ∈ (0.48 ∼ 0.84) fm.For small b ⊥ (b ⊥ ≤ 3a), we have further removed two points near the source and sink(t = 1 and t = t seq − 1) .
In Fig. 5, we give the lattice data and the fit results of the real and imaginary parts of C γ t 3 (t, t s )/C 2 (t s ) with different {P z , m π }, z and b ⊥ as examples.As shown in the figures, the fit results (colored bands) reproduce the original lattice data points at each t s , and the grey band corresponds to the extracted ground-state matrix element.
The fitting qualities of the bootstrap samples are illustrated in Fig. 6.The left panel shows the histogram distributions (normalized to 1) of χ 2 /d.o.f., and the right panel shows the cumulative distribution function (CDF) of the Q values for all correlated joint fits.Most of the χ 2 /d.o.f.spread out between 0.5 and 1.2, and the Q value is larger than 0.05 for most fits, which indicates that the ground-state fits are reasonable.
To demonstrate the stability of ground state fits, the comparisons between fit results of the bare matrix elements h0 γ t with different t min are shown in Fig. 7.The t min represents the minimum t s in the fits of the threepoint function, and the Q in the figure is the p-value of the joint fits.The two upper subdiagrams are at b ⊥ = 3a and z = 3a, for the real (left) and imaginary (right) part, respectively.The two lower panels are at b ⊥ = 5a and z = 5a.It shows that when t min increases, that is, when fitting with fewer data points, the fit results get larger uncertainties but remain consistent within the error range.

C. Renormalization
As mentioned in Sec.II B, we use the square root of the Wilson loop √ Z E and logarithmic divergence factor Z O to renormalize the bare quasi-TMD matrix elements.In practice, the signal to noise ratio of Z E (r, b ⊥ , a) defined in Eq.( 7) decreases fast with r and b ⊥ such that it is hardly available for large r or b ⊥ .To address this, we fit the effective energies of the Wilson loop, which gives the QCD static potentials, and then extrapolate them to large r and/or b ⊥ , as in Ref. [28].Numerical results for the Wilson loop are shown in the upper panel of Fig. 8.
The logarithmic divergence factor Z O can be extracted from the matching between lattice and perturbative re- sults in the region where both two theories work well.To preserve a good convergence of the perturbation theory before and after RG evolution, we choose the region where b ⊥ = 1a, and z = 0 or 1a, where both perturbation theory and lattice calculations work.The extracted Z O values at these two points are 1.0734(93) and 1.0509(92), respectively.Averaging these data points yields an aggregated result of Z O = 1.0622(87).
With Eq.( 8), it can be concluded that after dividing the bare matrix elements h0 Γ by √ Z E and Z O , the renormalized matrix elements should approximately be equal to the RG evolved perturbation results hMS Γ .The lower panel of Fig. 8 shows the consistency at points where we extract the Z O factor, which serves as a check of the numerical result for Z O .The points at z = 0 and z = 1 depicted in Fig. 8 demonstrate that, following renormalization with the extracted factor Z O = 1.0622(87), the renormalized matrix elements on the lattice show a consistent trend with the perturbative results obtained in MS scheme at short distance.To verify the accuracy of scale evolution from the intrinsic physical scale in lattice calculation to MS scale µ = 2GeV, we vary the origin of the evolution from 0.8µ 0 (shown as "RG.I") to 1.2µ 0 (shown as "RG.II"), which is sensitive to higher-order corrections.One can see that the perturbative results from different µ 0 are consistent with the lattice data in the matching region (b ⊥ = 1a, z = 0 or 1a).

D. L-dependence of subtracted quasi TMD matrix elements
For a well-defined quasi TMDPDF, the length of the Wilson link L should be large enough to ensure that the final results are independent of L. In the definition of the staple-shaped Wilson link in Eq.( 4), L should be extended to infinity.The L-dependence is included in the linear divergence of the Wilson link self-energy and pinchpole singularity, which can be subtracted by the square root of the Wilson loop.Therefore, the subtracted quasi TMD matrix elements will saturate to a plateau at large enough L to ensure that the link can extend outside the region of the parton, and exhibit independence of L.
Ref. [28] explored the L-dependence of renormalized TMD matrix elements and reached the conclusion that setting L = 6a for MILC12 is sufficient, aligning well with the ensemble we have employed.To verify this, we further examine the L-dependence of the subtracted quasi TMD matrix elements, illustrated in Fig. 9 and Fig. 10.
In Fig. 9, the three panels on the left show the real parts of the ratios defined in Eq.( 17 From the comparison, one can see that the fitted results for L = {6, 8, 10}a are consistent with each other.However, it should be noted that for large spatial separations, the signal becomes worse with increasing L. In Fig. 10, a comparison of the subtracted quasi TMDPDF as a function of λ with different L, also confirms this behavior.Therefore, to balance the statistical and systematic uncertainties, we adopt L = 6a as an optimal choice in our calculation. There is a phenomenological explanation for the earlier saturation of the Wilson line; however, it is essential to note that a definitive proof is currently unavailable.When calculating the TMDPDF of a proton, which reflects the quark correlation within the proton, a reasonable estimate for the correlation length could be the size of a proton, approximately 1 fm.Beyond this distance, quarks and gluons may escape the proton, potentially diminishing their impact.This suggests a typical saturation length of around 1 fm, with L = 6a (or more precisely, L = 8a) being in close proximity to this length scale.

E. Quasi TMDPDFs in the coordinate space and λ extrapolation
Combining the bare quasi TMDPDFs matrix elements with the corresponding Wilson loop and renormalization factor, we obtain numerical results for renormalized matrix elements at different λ = zP z .The two panels on the left in Fig. 11 exhibit the λ dependence of the renormalized matrix elements with b ⊥ = 2a and 5a with various P z and m π .It can be seen that as λ increases, the quasi TMDPDFs approach zero for both real and imaginary components.When b ⊥ is large, considerable uncertainties exist in the fully correlated datasets and non-zero central values can induce unphysical oscillations for a brute-force Fourier transformation.
To address this, we adopt a physics-inspired extrapolation at large λ: [32] hΓ,extra in which all parameters m 1,2 , n 1,2 and λ 0 depend on the transverse separation b ⊥ .The algebraic terms account for a power law behavior in the endpoint region, and the exponential term is motivated by the expectation that the correlation function has a finite correlation length (denoted as λ 0 ) at finite momentum.
It should be noticed that the λ extrapolation with Eq. ( 19) was was initially proposed for the one-dimensional parton distribution functions [32].For TMDPDFs, whether this form holds and whether the involved parameters depend on the transverse separation remain to be found out.In this work, we have performed the λ extrapolation in two ways: a joint extrapolation with the same parameters and an independent extrapolation for each b ⊥ .The results are collected in Tab.I, from which one can see that each separate fit result is consistent with the joint one.However, the joint fit approach gives a more stringent constraint for the large b ⊥ case, and the corresponding result for the extrapolated matrix element has a smaller error than the original lattice data, as shown in Fig. 12.To be conservative, we have adopted the independent fits for our results.
To perform the extrapolation, a reasonable range of λ is required to determine the parameters.In Fig. 13, we show several extrapolations of the quasi TMD matrix elements subtracted with different combinations of {m π , P z , b ⊥ }.The fit is performed in the region λ ≥ λ L with λ L being a truncation parameter.As one can see from the figures, the fitting bands agree with the original lattice data in the moderate λ region and have smooth tails at large λ.With that, the quasi TMDPDFs in momentum space do not have unphysical oscillations which often show up in the brute-force Fourier transformation of λ.As a conclusion, we choose λ L = 8a for the extrapolation and use λ L = 6a to estimate the systematic uncertainties from λ extrapolation.As shown in Fig. 11(c tice data in the moderate λ region, and give smoothlydecaying distributions at large λ.

F. TMDPDFs and physical extrapolation
After Fourier transforming the renormalized quasi TMDPDFs in coordinate space to momentum space, one can obtain the quasi TMDPDF fΓ (x, b ⊥ , ζ z , µ), shown as the dashdotted (magenta) line in Fig. 14.Combing the updated results of Collins-Soper kernel [27] and intrinsic soft function [40] calculated on the same ensembles, we obtain the TMDPDFs through the matching formula in Eq. (10).Fig. 14 shows an example of matched TMD-PDF from the NLO [22,42] and NNLO [43,44] kernel with RG running to scale µ = √ ζ = 2 GeV.One can see that the results agree within uncertainties, except for the endpoint region.
It should be noted that the TMDPDFs can be obtained after employing the factorization formula, while residual dependence on the lattice inputs (such as P z and m π ) may still reside in the obtained results.To diminish the dependence, we extrapolate the above results to the physical m π value (135 MeV) and infinite momentum using the following ansatz: where the d 0 term characterizes the pion mass dependence, and d 1 accounts for the momentum-dependent discretization error.An interesting analysis has recently derived the chiral logarithms for quasi-PDFs [51].In this approach, one first performs an operator product expansion and constructs the moments of quasi PDFs  I. Results for the parameters n1, n2 and λ in Eq.( 19) from separate fits and a joint fit, and the χ 2 /d.o.f. of each fit, taking the case of mπ = 220 MeV, Γ = γ t and P z = 1.72GeV as an example.The results of different fit methods are consistent with each other.
! ℎ with derivative operators.Using chiral perturbation theory, one can find hadron-level operators which have the same symmetry.Calculating the one-loop perturbative contributions and summing these contributions for moment operators will lead to the chiral logarithms in quasi-PDFs.However, this analysis has not been generalized to TMDPDFs at present.In addition, since TMDPDFs are three-dimensional distributions, expanding the nonlocal operators will require the derivatives not only in the z direction, but also in the transverse direction.This will likely introduce additional complexities.
Taking d ′ 0 = d ′ 1 = 0 as our default case, we obtain the physical TMDPDF from a joint fit of f (m π , P z ) shown in Fig. 15.In order to explore the systematic bias in the extrapolation ansatz, we also estimate the m π dependence by using the d ′ 0 term (with d 0 = 0) and examine the 1/P z contributions by adding the d ′ 1 term, shown as the following two different strategies: • Default form: • Strategy I: Specifically, we introduce the chiral log term as well as the linear 1/P z term into the fit formula, and consider the between them and the default form as the systematic uncertainty.A comparison of them is shown in the upper panel in Fig. 16 from which one can see that the results agree with each other.Furthermore, we have examined two different sequences of extrapolation: First, extrapolation to the physical m π followed by the P z → ∞ extrapolation, and second, the P z → ∞ extrapolation followed by the extrapolation to the physical m π .A comparison of the results is presented in the lower panel of Fig. 16.It is evident that the results are in good agreement with each other.
G. γ t − γ z difference In quasi TMDPDFs, both Lorentz structures Γ = γ t and γ z can project onto γ + in the large momentum limit.An interesting observation in Ref. [52] is that the γ t case has fewer operator mixing effects.Therefore, we use the Γ = γ t in Eq.( 1) to extract the TMDPDFs.For the γ z case, in addition to the sizable operator mixing effects, it is anticipated that its deviations from the γ t results may come from power corrections from the operator product expansion of quasi correlators.These corrections are of order O(1/(P z ) 2 ) with opposite signs at large P z .In order to analyze the impact caused by different structures, we plot the ratio |f γ t − f γ z |/(f γ t + f γ z ) with different P z in Fig. 17.One can see that with increasing nucleon momentum, the ratio becomes smaller, except for the endpoint region with large uncertainties.The differences will be incorporated as a systematic uncertainty.

H. Estimation of systematic uncertainties
As mentioned above, this work considers the systematic uncertainties from different sources, including from • different fitting ranges in the λ-extrapolation (Sec.III E); • different strategies used for the chiral and P z extrapolation (Sec.III F); • the difference between γ t and γ z results (Sec.III G).
In addition, we also consider the error propagating from the intrinsic soft function [40] and Collins-Soper kernel [27], which are calculated on the same configurations.
The comparison of the statistic as well as each systematic uncertainties are shown in Fig. 18.

IV. FINAL RESULTS FOR TMDPDFS
Combing all the known uncertainties indicated in Fig. 18, we obtain numerical results of unpolarized nucleon's isovector TMDPDFs from our lattice simulation.In Fig. 19, a bump in the x distribution can be observed, representing the highest probability for the partons' longitudinal momentum distributions.It is worth noting that the peak positions do not align precisely between the lattice results and the respective phenomenological results.Further investigations are required to elucidate this phenomenon.
The two shaded bands at the endpoint regions (x < 0.2 and x > 0.8) in subplots of Fig. 19 indicate that LaMET predictions are not reliable there, which is estimated from the power correction terms Λ 2 QCD /(xP z ) 2 and threshold logarithms ln((1 − x)P z ) [53,54] .Besides, since only one lattice spacing is used, discretization uncertainties are not properly handled at this stage.Especially the b ⊥ ∼ a case might suffer sizable discretization effects, which can be improved by a more detailed analysis in future.GeV, extrapolated to physical pion mass 135 MeV and infinite momentum limit P z → ∞, compared with ART23 [7], BHLSVZ22 [8], MAPTMD22 [9], SV19 [10]and PV17 [11] global fits.The colored bands denote our results with both statistical and systematic uncertainties, the shaded grey regions imply the endpoint regions where LaMET predictions are not reliable.

V. SUMMARY AND PROSPECT
In summary, we have performed the first calculation of Transverse Momentum-Dependent Parton Distribution Functions (TMDPDFs) within a nucleon using the LaMET expansion of lattice data.State-of-the-art techniques in renormalization and extrapolation on the lattice have been employed, including the consideration of the perturbative kernel up to NNLO with RG resummation.We investigate the dependence on pion mass and hadron momentum, incorporating both statistical and systematic errors to provide a robust characterization of the inner structure of nucleons through the parton distributions.
While the current findings are promising, further improvements are required.Firstly, a comprehensive anal-ysis of multiple ensembles with varying lattice spacings, pion masses, and volumes is essential.This approach will only reassess all uncertainties addressed in study but also systematically explore additional significant factors.Secondly, our results exhibit notable theoretical uncertainties in the endpoint regions, which could stem from power corrections or threshold logarithms.Thirdly, we expect a decay of the TMDPDF with increasing b ⊥ which is not (yet) visible in our data.Going to larger b ⊥ is thus imperative.Moreover, the presence of chiral logarithms in the analysis of quasi PDFs (as demonstrated in Ref. [51]) highlights the potential for similar logarithms in our case, warranting dedicated theoretical investigation.

FIG. 2 .
FIG. 2. Comparison of NLO and NNLO matching kernels at fixed order with µ = 2GeV (dashed lines) and running from Collins-Soper scale µ0 = 2xP z to MS scale µ = 2GeV (solid lines), P z is chosen as 2.15GeV.The abrupt surge in behavior observed in NNLO+RGR is attributed to the presence of the Landau pole.

) ]FIG. 5 .FIG. 6 .
FIG. 5. Ratios of C γ t3 (t, ts)/C2 (ts) (data points), as functions of t and ts, with various combination of {P z , mπ} and different z and b ⊥ .In this figure, the colored bands correspond to the fitted results, and the grey bands are the ground-state contribution.

FIG. 7 .
FIG.7.Stability plots to compare fit results of bare matrix elements when varying the minimum ts in the fits of C γ t 3 (t, ts)/C2 (ts).The two upper subdiagrams denote the real and imaginary results at (b ⊥ , z) = (3, 3)a and the lower ones are at (5, 5)a.The Q in each subplots represents the p-value of the fits, all the fits above have good qualities with p-values greater than 0.1, which are denoted with the red dashed lines.
FIG. 9.The real parts of ratios defined in Eq.(17) divided by the Wilson loop (three left panels) and the L dependence of the ground state results after fitting (right panels) at {z, b ⊥ } = {1, 1}, {3, 3} and {5, 5}a, with P z = 1.72 GeV, mπ = 220 MeV and Γ = γ t .Since the real parts for {z, b ⊥ } = {5, 5}a are close to zero, we also give the imaginary parts in this example.The labels are similar to the corresponding ones in Fig.5.

𝑃 ! = 1 FIG. 15 .
FIG. 14. Quasi TMDPDF (dashed-red line) with nucleon boosted momentum P z = 1.72 GeV and matched TMDPDF from NLO kernel (dashed-black line, no error) and NNLO kernel (solid-black line) with RG running to scale µ = √ ζ = 2 GeV at b ⊥ = 3a.Only statistical errors are included in the bands.Deviations between NLO and NNLO results at x < 0.2 indicate that the perturbative matching fails in the small-x region.
FIG.16.Upper panel: Default extrapolation form in Eq.(21), considering the systematic uncertainties from chiral extrapolation in Eq.(22) (strategy I) and from the large-momentum extrapolation form in Eq.(23) (strategy II).Lower panel: Comparison of the results from joint fit and chained fits with different order.

Fig. 19
Fig. 19 shows xf (x, b ⊥ , µ, ζ) at renormalization and rapidity scales µ = √ ζ = 2 GeV as a function of x, together with the phenomenological results from global analyses [7-11].From the comparison one can see that, our results are in qualitative agreement with phenomenological results and share a similar behavior in b ⊥ space: the central values slowly decrease and uncertainties are gradually increasing with the increase of b ⊥ .

Fig. 20 6 ART23FIG. 19 .
Fig.20shows the results for xf (x, b ⊥ , µ, ζ) with x = 0.3, 0.5, 0.7 as a function of spatial separation b ⊥ .The spatial distributions reflect correlations between the partons at transverse interval b ⊥ inside a highly boosted nucleon, and will reveal aspects of nucleon structure.It could be conjectured that the distributions vanish when b ⊥ is larger than the nucleon radius, however, neither the present day lattice results nor the phenomenological results are precise enough to draw such a conclusion.