Transverse-momentum-dependent pion structures from lattice QCD:
Collins-Soper kernel, soft factor, TMDWF, and TMDPDF

Abstract

We present the first lattice quantum chromodynamics (QCD) calculation of the pion valence-quark transverse-momentum-dependent parton distribution function (TMDPDF) within the framework of large-momentum effective theory (LaMET). Using correlators fixed in the Coulomb gauge (CG), we computed the quasi-TMD beam function for a pion with a mass of 300 MeV, a fine lattice spacing of $a=0.06$ fm and multiple large momenta up to 3 GeV. The intrinsic soft functions in the CG approach are extracted from form factors with large momentum transfer, and as a byproduct, we also obtain the corresponding Collins-Soper (CS) kernel. Our determinations of both the soft function and the CS kernel agree with perturbation theory at small transverse separations ($b_{\perp}$) between the quarks. At larger $b_{\perp}$, the CS kernel remains consistent with recent results obtained using both CG and GI TMD correlators in the literature. By combining next-to-leading logarithmic (NLL) factorization of the quasi-TMD beam function and the soft function, we obtain $x$-dependent pion valence-quark TMDPDF for transverse separations $b_{\perp}\gtrsim 1$ fm. Interestingly, we find that the $b_{\perp}$ dependence of the phenomenological parameterizations of TMDPDF for moderate values of $x$ are in reasonable agreement with our QCD determinations. In addition, we present results for the transverse-momentum-dependent wave function (TMDWF) for a heavier pion with 670 MeV mass.

I Introduction

In high-energy scatterings, transverse-momentum-dependent parton distribution functions (TMDPDFs) provide a fundamental description of the transverse momentum and polarization degrees of freedom of quarks and gluons within hadrons [1]. These distributions play a crucial role in understanding the intricate dynamics of quark-gluon interactions and the phenomenon of color confinement. Knowledge of TMDPDFs is essential for predicting observables in multi-scale, non-inclusive high-energy processes, such as semi-inclusive deep-inelastic scattering (SIDIS) and Drell-Yan (DY), based on QCD factorization. Experiments at facilities including the Jefferson Lab 12 GeV Program [2], the Electron-Ion Collider (EIC) [3, 4], and the Large Hadron Collider (LHC) [5] rely heavily on accurate knowledge of the TMDPDFs.

Significant progress has been made in the phenomenological parameterizations of TMDPDFs, particularly for the nucleon, through global analyses of experimental data [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. While these studies have significantly improved our understanding of the transverse-momentum structure of quarks and gluons within the nucleon, they remain in their early stages due to the limited availability of experimental data sensitive to the non-perturbative region. Compared to nucleons, much less is known about the transverse momentum structure of the pion [27, 28]. As the lightest QCD bound state, the pion plays a fundamental role in hadron structure and non-perturbative QCD. Its TMDPDFs provide crucial insights into the internal dynamics of quarks and gluons in a strongly interacting relativistic system, with direct implications for hadronization processes and the underlying mechanisms of non-perturbative QCD. However, the scarcity of experimental data on pion TMDPDFs, combined with the inherent nature of the parameterization dependence of global analyses, leaves significant gaps in our non-perturbative first-principles understanding. Lattice QCD calculations provide a natural approach to gain insight into non-perturbative TMD structures of pion starting from first-principles QCD.

Early lattice QCD studies primarily focused on computing the moments of TMDPDFs [29, 30, 31, 32]. More recently, large-momentum effective theory (LaMET) [33, 34, 35, 36] has provided a framework for directly calculating the $x$-dependence of parton distributions [37], including TMDPDFs [38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55]. In the large-momentum limit, quasi-TMDs, defined as equal-time correlators in highly boosted hadron states, can be related to their light-cone counterparts through perturbative factorization, making first-principles calculations of TMDPDFs possible within lattice QCD.

This framework has driven substantial progress in recent years. One of its key achievements is the determination of the Collins-Soper (CS) kernel, which governs the scale dependence of TMDPDFs [56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66]. Another crucial development is the proposal to extract the intrinsic soft function [43, 52] - which accounts for soft gluon radiation and plays a vital role in TMD factorization - from meson form factors with large momentum transfer [57, 58, 63]. This step completes the factorization framework that connects quasi-TMDs to light-cone TMDs in the large-momentum limit. These advances have significantly enhanced first-principles lattice calculations of TMDPDFs. Recent achievements include studies of the pion TMD wave function (TMDWF) [67], the extraction of the unpolarized nucleon TMDPDF [68], and investigations of the Boer-Mulders functions of the pion and nucleon [69, 70].

Despite these advancements, a major challenge in lattice calculations of gauge-invariant (GI) quasi-TMDPDFs arises from the structure of the quark-bilinear operators, which are separated in both longitudinal and transverse directions and connected by a staple-shaped Wilson link to maintain gauge invariance. The presence of Wilson lines introduces linear divergences, requiring careful renormalization [71, 72, 42, 73, 74, 64, 75, 76], and leads to an exponential suppression of the signal-to-noise ratio (SNR) as the quark field separation increases [66]. This suppression presents a significant obstacle to precise lattice determinations of TMDPDFs, particularly in the low transverse momentum region or at large spatial separations, where high precision is crucial. Furthermore, controlling the discretization effects and power corrections at finite momentum remains a challenge in ensuring reliable extrapolations to the continuum and physical limits.

To address this issue, a new approach based on LaMET has been proposed [77, 78] to extract TMDPDFs from Coulomb gauge (CG) fixed quark correlators. Unlike traditional GI correlators, CG correlators do not require Wilson lines. Despite this difference, they fall within the same universality class [79, 35] as GI correlators with staple-shaped Wilson lines, as both reduce to light-cone TMD correlators in the infinite boost limit. Consequently, CG quasi-TMDPDFs can be matched to light-cone TMDPDFs through a perturbative factorization [78]. The absence of Wilson lines simplifies renormalization and significantly enhances the SNR of boosted correlators, particularly at large transverse separations.

In this study, we present the first lattice QCD calculation of the unpolarized valence-quark pion TMDPDF using the CG method. We compute the quasi-TMDPDF and quasi-TMDWF matrix elements of the pion in CG. These matrix elements enable the extraction of the CS kernel, the intrinsic soft function, and ultimately the unpolarized valence TMDPDF of the pion. Furthermore, we improve the perturbative accuracy of the matching procedure, particularly in the extraction of the intrinsic soft function, by resumming logarithmic terms in the Sudakov kernel, an aspect that has not been accounted for in previous lattice QCD studies. Finally, the comparison of our results with global fits of the experimental data shows encouraging consistency, demonstrating the bright potential of the CG approach for TMD physics from lattice QCD.

The paper is organized as follows: we first introduce the theoretical framework for calculating the TMDPDF from the CG quasi-TMDPDFs in Sec. II, which also includes the CS kernel and the intrinsic soft factor; the lattice setup is presented in Sec. III; then we present the detailed analysis of the quasi-TMD in Sec. IV; the CS kernel, the intrinsic soft factor, the full pion TMDWF and TMDPDF are analyzed in Sec. V to be compared with phenomenological results; finally, we conclude in Sec. VI.

II Theoretical framework

II.1 Light-cone TMDPDF from CG quasi-TMDPDF

As proposed in Ref. [78], the light-cone TMDPDF can be derived from the CG quasi-TMD beam function defined as,

$$\tilde{f}_{\Gamma}(x,b_{\perp},P^{z};\mu)=P^{z}\int\frac{dz}{2\pi}e^{iz(xP^{z})}\tilde{h}_{\Gamma}(z,b_{\perp},P^{z};\mu)~{},$$ (1)

where the matrix elements are given by,

$$\tilde{h}_{\Gamma}(z,b_{\perp},P^{z};\mu)=\frac{1}{2P^{t}}\bra{P}\left.\bar{q}(z,b_{\perp})\Gamma q(0)\right|_{\vec{\nabla}\cdot\vec{A}=0}\ket{P}~{}.$$ (2)

In this work, we choose $\Gamma=\gamma^{t}\gamma^{5}$ for quasi-TMD beam function. The parameter $\mu$ represents the $\overline{\text{MS}}$ renormalization scale, and the space-like separation between the quark and the antiquark is denoted as $\vec{b}\equiv(z,\vec{b}_{\perp})$. The hadron state carries momentum $P=(P^{t},0,0,P^{z})$ and satisfies the normalization condition $\langle P|P\rangle=2P^{t}\delta(0)$. When the hadron moves with a large momentum $P^{z}$, the quasi-TMD beam function factorizes as,

$$\begin{split}\tilde{f}_{\Gamma}(x,b_{\perp},P^{z};\mu)=H_{f}&\left(x,P^{z};\mu\right)B(x,b_{\perp},xP^{+};\mu,\nu)\\ &\times S^{0}_{C}(b_{\perp};\mu,\nu)~{},\end{split}$$ (3)

where $B(x,b_{\perp},xP^{+};\mu,\nu)$ is the light-cone beam function with zero-bin subtraction, following the procedure outlined in Refs. [80, 81]. The term $S^{0}_{C}(b_{\perp};\mu,\nu)$ represents the zero-bin contribution, and the operator definitions of both functions are detailed in Ref. [78]. The parameters $\mu$ and $\nu$ correspond to the renormalization scales associated with the ultraviolet (UV) and rapidity divergences, respectively. The function $H_{f}\left(x,P^{z};\mu\right)=|C_{\rm TMD}(xP^{z};\mu)|^{2}$ represents the hard kernel that matches the QCD quark field operator to Soft-Collinear Effective Theory (SCET) [45].

The light-cone TMDPDF can be defined from the light-cone beam function $B$ and the soft function $S$,

$$f(x,b_{\perp};\mu,\zeta)=B(x,b_{\perp},xP^{+};\mu,\nu)S(b_{\perp},y_{n};\mu,\nu)~{},$$ (4)

where the dependence on the rapidity renormalization scale $\nu$ cancels out on the right-hand side, leaving only a dependence on the Collins-Soper scale $\zeta\equiv 2(xP^{+})^{2}e^{-2y_{n}}$. Combining Eqs. 3 and 4, the factorization formula connecting the quasi-TMDPDF and the light-cone TMDPDF is given by [40, 41, 43, 44, 78]

$$\begin{aligned} &\sqrt{S_{I}\left(b_{\perp},y_{n};\mu\right)}\cdot\tilde{f}_{\Gamma}(x,b_{\perp},P^{z};\mu)\\ &\quad\quad\quad\quad\quad=f(x,b_{\perp};\mu,\zeta)H_{f}\left(x,P^{z};\mu\right)+\mathrm{p.c.}~{},\end{aligned}$$ (5)

where p.c. indicates the power corrections, and the intrinsic soft function $S_{I}$ is defined as

$$S_{I}\left(b_{\perp},y_{n};\mu\right)\equiv\left(\frac{S(b_{\perp},y_{n};\mu,\nu)}{S^{0}_{C}(b_{\perp};\mu,\nu)}\right)^{2}~{}\,.$$ (6)

Here $S^{0}_{C}$ is the CG quasi-soft function, as defined in Ref. [78], which exactly cancels the rapidity divergences in $S$. The above factorization formula is of the same form as that for the GI quasi-TMDPDFs [40, 41, 43, 44], while the method to calculate $S_{I}$ was first proposed in Ref. [43], which enables a complete determination of the TMDs from the lattice.

The intrinsic soft function satisfies the evolution equation

$$\frac{d}{dy_{n}}\ln S_{I}\left(b_{\perp},y_{n};\mu\right)=-2\gamma^{\overline{\rm MS}}(b_{\perp};\mu)~{},$$ (7)

where $\gamma^{\overline{\rm MS}}(b_{\perp};\mu)$ is the Collins-Soper kernel. Therefore, the intrinsic soft function at any $y_{n}$ satisfies

$$\begin{aligned} &S_{I}\left(b_{\perp},y_{n};\mu\right)\\ &\quad=S_{I}\left(b_{\perp},y_{n}=0;\mu\right)\cdot\exp(-2y_{n}\cdot\gamma^{\overline{\rm MS}}(b_{\perp};\mu))~{}.\end{aligned}$$ (8)

One can redefine the intrinsic soft function as $S_{I}\left(b_{\perp};\mu\right)\equiv S_{I}\left(b_{\perp},y_{n}=0;\mu\right)$, allowing us to rewrite Eq. (5) in a more explicit form,

$$\begin{aligned} &\sqrt{S_{I}\left(b_{\perp};\mu\right)}\cdot\tilde{f}_{\Gamma}\left(x,b_{\perp},P^{z};\mu\right)=f\left(x,b_{\perp};\mu,\zeta\right)\\ &\quad\quad\times H_{f}\left(x,P^{z};\mu\right)\exp\left[\frac{1}{2}\ln\frac{\left(2xP^{z}\right)^{2}}{\zeta}\gamma^{\overline{\mathrm{MS}}}\left(b_{\perp};\mu\right)\right]\\ &\quad\quad+\mathcal{O}\left(\frac{\Lambda_{\mathrm{QCD}}}{xP^{z}},\frac{1}{b_{\perp}\left(xP^{z}\right)}\right)~{}.\end{aligned}$$ (9)

This factorization framework provides a systematic approach to relating quasi-TMDPDFs extracted from lattice QCD to their light-cone counterparts, ensuring a well-defined separation between perturbative and non-perturbative contributions. However, power corrections of order $\mathcal{O}\left(\frac{\Lambda_{\mathrm{QCD}}}{xP^{z}},\frac{1}{b_{\perp}xP^{z}}\right)$ [82] introduce systematic uncertainties, particularly in the small-$x$ region. Despite these limitations, the framework offers valuable first-principles constraints in the large-$x$ regime, where power corrections are expected to be relatively small.

II.2 Quasi-TMDWF and CS kernel

Similar to the quasi-TMD beam function defined in Eqs. (1) and (2), the light-cone TMDWF can be derived from the CG quasi-TMDWF, which is defined as

$$\tilde{\phi}_{\Gamma}(x,b_{\perp},P^{z};\mu)=P^{z}\int\frac{dz}{2\pi}e^{iz(xP^{z})}\tilde{\varphi}_{\Gamma}(z,b_{\perp},P^{z};\mu)~{},$$ (10)

with the quasi-TMDWF matrix elements given by

$$\begin{aligned} &\tilde{\varphi}_{\Gamma}(z,b_{\perp},P^{z};\mu)=\frac{e^{-iz(xP^{z})/2}}{2P_{t}}\\ &\quad\times\bra{\Omega}\left.\bar{q}(z/2,b_{\perp}/2)\Gamma q(-z/2,-b_{\perp}/2)\right|_{\vec{\nabla}\cdot\vec{A}=0}\ket{P}~{},\end{aligned}$$ (11)

which corresponds to a pion-to-vacuum matrix element, the $\ket{\Omega}$ represents the QCD vacuum state. In this work, we choose $\Gamma=\gamma^{z}\gamma^{5}$ for quasi-TMDWFs with nonzero momenta, while for zero-momentum quasi-TMDWFs, we set $\Gamma=\gamma^{t}\gamma^{5}$.

In the large-momentum limit, the quasi-TMDWF can be matched to the light-cone TMDWF via the factorization formula

$$\begin{aligned} &\sqrt{S_{I}\left(b_{\perp},y_{n};\mu\right)}\cdot\tilde{\phi}_{\Gamma}\left(x,b_{\perp},P^{z};\mu\right)\\ &\quad\quad=\phi(x,b_{\perp},y_{n};\mu,\zeta,\bar{\zeta})H_{\phi}\left(x,\bar{x},P^{z};\mu\right)+\mathrm{p.c.}~{},\end{aligned}$$ (12)

with $\bar{x}=1-x$, $\bar{\zeta}\equiv 2(\bar{x}P^{+})^{2}e^{-2y_{n}}$, and the hard kernel

$$H_{\phi}(x,\bar{x},P^{z};\mu)=C_{\rm TMD}(xP^{z};\mu)\cdot C_{\rm TMD}(\bar{x}P^{z};\mu)~{},$$ (13)

calculated up to NLO in Ref. [78].

An essential component of quasi-TMD factorization is the CS kernel $\gamma^{\overline{\rm MS}}(b_{\perp};\mu)$, which governs the rapidity evolution of TMDs and is crucial to achieve consistent matching between quasi-TMDPDFs and light-cone TMDPDFs. The CS kernel can be extracted from the ratios of quasi-TMD matrix elements computed at different hadron momenta [38, 40, 43]. In this work, we employ the ratio of CG quasi-TMDWFs [43, 48],

$$\begin{aligned} &\gamma^{\overline{\rm MS}}(b_{\perp};\mu)\\ &\quad=\gamma^{\overline{\rm MS}}(b_{\perp},P_{1},P_{2};\mu)+\text{p.c.}\\ &\quad=\frac{1}{\ln\left(P_{2}/P_{1}\right)}\ln\frac{H_{\phi}\left(x,\bar{x},P_{1};\mu\right)\tilde{\phi}_{\gamma^{z}\gamma^{5}}\left(x,b_{\perp},P_{2};\mu\right)}{H_{\phi}\left(x,\bar{x},P_{2};\mu\right)\tilde{\phi}_{\gamma^{z}\gamma^{5}}\left(x,b_{\perp},P_{1};\mu\right)}\\ &\quad\quad+\mathcal{O}\left(\frac{\Lambda_{\mathrm{QCD}}}{xP^{z}},\frac{1}{b_{\perp}\left(xP^{z}\right)},\frac{\Lambda_{\mathrm{QCD}}}{\bar{x}P^{z}},\frac{1}{b_{\perp}\left(\bar{x}P^{z}\right)}\right)~{}.\end{aligned}$$ (14)

This method offers a key advantage over direct extractions from quasi-TMDPDFs, as quasi-TMDWFs naturally peak around $x=0.5$, a region where power corrections are significantly suppressed. On the other hand, the quasi-TMDPDFs usually decreases quickly at moderate to large $x$, making the extraction of the CS kernel more susceptible to power corrections.

II.3 Intrinsic soft function

Another essential component of the factorization of quasi-TMDs is the intrinsic soft function, which accounts for soft gluon radiation in the process. Its calculation is particularly challenging for two reasons. First, as a non-perturbative quantity, it cannot be computed using standard perturbative techniques. Second, because the soft function involves correlations along two light-cone directions, it cannot be directly simulated on the lattice, even with a large-momentum boost.

A practical approach to determining the intrinsic soft function, as proposed in Refs. [43], is to use the TMD factorization of a pseudoscalar light-meson form factor, which is defined as

$$\begin{aligned} &F\left(b_{\perp},P_{1},P_{2},\Gamma,\Gamma^{\prime},\Gamma^{*}\right)\\ &\quad\equiv\frac{\left\langle P_{2}\right|\bar{q}\left(b_{\perp}\right)\Gamma q\left(b_{\perp}\right)\bar{q}(0)\Gamma^{\prime}q(0)\left|P_{1}\right\rangle}{\langle 0|\bar{q}(0)\Gamma^{*}q(0)\left|P_{1}\right\rangle\left\langle P_{2}\right|\bar{q}(0)\Gamma^{*}q(0)|0\rangle}~{},\end{aligned}$$ (15)

where the numerator represents a three-point function, with $\ket{P_{1}}$ and $\ket{P_{2}}$ denoting pion states separated by a time interval $t_{\rm sep}$ with momenta $P_{1}$ and $P_{2}$. The two current operators, $\bar{q}\left(b_{\perp}\right)\Gamma q\left(b_{\perp}\right)$ and $\bar{q}(0)\Gamma^{\prime}q(0)$, are inserted at the same time slice $\tau$ between these pion states. To extract the leading-twist contribution, one can choose the Dirac matrices as $\Gamma=\Gamma^{\prime}\in\{I,\gamma_{5},\gamma_{\perp},\gamma_{5}\gamma_{\perp}\}$ [52, 63]. The denominator consists of pion decay constants used to normalize the form factor, with the conventional choice of $\Gamma^{*}=\gamma^{t}\gamma^{5}$. The pion states are taken as $P_{1}=(P^{t},0,0,P^{z})$ and $P_{2}=(P^{t},0,0,-P^{z})$, moving back-to-back to achieve a large momentum transfer $Q^{2}=(2P^{z})^{2}$. In this kinematic regime, the form factor can be factorized in terms of the TMDWF $\phi(x,b_{\perp},y_{n};\mu,\zeta)$ [43, 48, 52] as

$$\begin{aligned} &F(b_{\perp},P^{z})=\int dx_{1}dx_{2}H_{F}(x_{1},x_{2},P^{z};\mu)\\ &\quad\quad\times\phi^{\dagger}(x_{1},b_{\perp},y_{n};\mu,\zeta_{1},\bar{\zeta}_{1})\phi(x_{2},b_{\perp},-y_{n};\mu,\zeta_{2},\bar{\zeta}_{2})~{}.\end{aligned}$$ (16)

where $H_{F}$ is the hard kernel encapsulating the short-distance dynamics of the process. The $y_{n}$-dependence cancels between $\phi$ and $\phi^{\dagger}$, and the hard kernel can be expressed as the product of two Sudakov kernels [52]:

$$\begin{aligned} &H_{F}(x_{1},x_{2},P^{z};\mu)\\ &\quad=C_{\rm Sud}(x_{1},x_{2},P^{z};\mu)\cdot C_{\rm Sud}(\bar{x}_{1},\bar{x}_{2},P^{z};\mu)~{},\end{aligned}$$ (17)

where $C_{\rm Sud}$ represents the Sudakov kernel, which has been computed at one-loop order in the literature [83, 52]. The renormalization group (RG)-resummed results for $C_{\rm Sud}$ at next-to-leading logarithmic (NLL) accuracy are provided in App. A.

By combining the factorization of the quasi-TMDWF in Eq. (12) with the form factor factorization in Eq. (16), the intrinsic soft function at $y_{n}=0$ can be extracted as

$$S_{I}(b_{\perp};\mu)=\frac{F(b_{\perp},P^{z})}{\int dx_{1}dx_{2}H_{F}(x_{1},x_{2},P^{z};\mu)\tilde{\Phi}^{\dagger}(x_{1})\tilde{\Phi}(x_{2})}~{},$$ (18)

where the reduced quasi-TMDWF $\tilde{\Phi}(x)$, is defined as

$$\tilde{\Phi}(x)\equiv\frac{\tilde{\phi}_{\Gamma}\left(x,b_{\perp},P^{z};\mu\right)}{H_{\phi}\left(x,\bar{x},P^{z};\mu\right)}~{}.$$ (19)

Using different renormalization schemes for quasi-TMDWF, one can get different intrinsic soft functions, which can be perturbatively converted to each other at small $b_{\perp}$. However, the scheme dependence will eventually cancel between the renormalized quasi-beam and intrinsic soft functions. The details of scheme conversion can be found in App. B.

III Lattice setup

We perform a numerical lattice QCD calculation using a $2+1$-flavor gauge ensemble generated by the HotQCD Collaboration [84]. The ensemble employs the Highly Improved Staggered Quark (HISQ) action [85] with a lattice spacing of $a=0.06$ fm and a volume of $L_{s}^{3}\times L_{t}=48^{3}\times 64$. The valence sector is treated using tadpole-improved clover Wilson fermions on a hypercubic (HYP) smeared [86] gauge background. The clover coefficient is set to $c_{\rm sw}=u_{0}^{-3/4}$, where $u_{0}$ is the average plaquette after HYP smearing. For this ensemble, we use $c_{\rm sw}=1.0336$ in both the time and spatial directions. The valence quark masses are tuned to yield a valence pion mass of 300 MeV, with a corresponding hopping parameter of $\kappa\approx 0.12623$.

The key requirement for the factorization of quasi-TMDs is a sufficiently large hadron momentum. To achieve a higher boost factor, we take advantage of the three-dimensional rotational symmetry of the CG approach and adopt off-axis momenta for the quasi-TMD, choosing momentum directions along $\vec{n}=(n^{x},n^{y},0)$. The hadron momenta on the lattice are given by $P^{z}=\frac{2\pi|\vec{n}|}{L_{s}a}$. In this study, we consider $n^{x}=n^{y}\in\{3,4,5\}$, which allows us to reach a maximum hadron momentum of $P^{z}=3.04$ GeV for the quasi-TMD calculations. To optimize the signal-to-noise ratio and enhance overlap with large-momentum hadron ground states, we employ boosted Gaussian smearing [87], using the same setup as in Ref. [77]. To extract the ground-state contribution, we compute quasi-TMD three-point functions for multiple source-sink separations, choosing $t_{\rm sep}/a=6,8,10,12$. Calculations are performed on 553 gauge configurations and we apply the All-Mode Averaging (AMA) technique [88] to further improve the signal. The stopping criteria for the sloppy and exact inversions are set to $10^{-4}$ and $10^{-10}$, respectively, aligning with the settings in Ref. [89]. For the smaller separation $t_{\rm sep}/a=6$, we compute 1 exact source and 32 sloppy sources per configuration, while for the larger separations $t_{\rm sep}/a=8,10,12$, we compute 4 exact sources and 128 sloppy sources.

In this work, we do not compute the form factors directly, but instead incorporate the form factor results from Ref. [63], which are computed on MILC ensembles with $2+1+1$-flavor with HISQ action [85]. Since those form factor calculations were performed using a different valence pion mass of 670 MeV, we compute the quasi-TMDWF with the same valence pion mass to ensure consistency in the extraction of the intrinsic soft function, which is a vacuum property and should be insensitive to the valence quark masses. To maintain this consistency, we adopt the conventional momentum setting $\vec{n}=(0,0,n^{z})$ and compute for $n^{z}\in\{8,9,10\}$, corresponding to a maximum hadron momentum of $P^{z}=4.30$ GeV. This calculation is performed on the same 553 gauge configurations from HotQCD Collaboration as mentioned above, employing 1 exact source and 32 sloppy sources per configuration.

IV quasi-TMDPDF matrix elements

IV.1 Two-point function and dispersion relation

Figure 1: Upper panel: effective mass plot of $C^{\rm SS}_{\text{2pt}}$ with various momenta. The dashed lines are the ground state energies of pion calculated using the dispersion relation $E_{0}=\sqrt{m_{0}^{2}+(P^{z})^{2}}$, where the valence pion mass is $m_{0}=300$ MeV. Lower panel: the ground-state energies $E_{0}$ extracted from two-state fit of the two-point functions. The red line represents the exact dispersion relations with $m_{0}=300$ MeV.

To extract the ground-state quasi-TMD matrix elements, it is essential to first determine the energy spectrum associated with the pion interpolator $\pi^{\dagger}$. The two-point correlator is defined as

$$C^{ss^{\prime}}_{\text{2pt}}\left(t_{\rm sep},P^{z}\right)=\langle\pi_{s^{\prime}}(P^{z},t_{\rm sep})\pi^{\dagger}_{s}(P^{z},0)\rangle~{},$$ (20)

where $\pi^{\dagger}_{s}$ and $\pi_{s^{\prime}}$ represent the pion source and sink defined by,

$$\begin{aligned} &\pi_{s}(\vec{x},t_{\rm sep})=\overline{d}_{s}(\vec{x},t_{\rm sep})\gamma_{5}u_{s}(\vec{x},t_{\rm sep}),\\ &\pi_{s}(\vec{P},t_{\rm sep})=\sum_{\vec{x}}\pi_{s}(\vec{x},t_{\rm sep})e^{-i\vec{P}\cdot\vec{x}}.\end{aligned}$$ (21)

In this work, we use a Gaussian-smeared source and sink ($s=s^{\prime}=S$), following the setup in Ref. [90]. The superscript $SS$ will be removed in the following text for simplicity.

By inserting a complete set of states, the two-point correlator can be expressed as a sum over energy eigenstates,

$$C_{\text{2pt}}\left(t_{\rm sep}\right)=\sum_{n=0}^{N_{\mathrm{s}}-1}\frac{\left|z_{n}\right|^{2}}{2E_{n}}\left(e^{-E_{n}t_{\rm sep}}+e^{-E_{n}\left(L_{t}-t_{\rm sep}\right)}\right)~{}.$$ (22)

Here, the overlap amplitude $z_{n}=\langle n|\pi^{\dagger}|\Omega\rangle$ quantifies the projection of the pion interpolator onto the $n$th energy eigenstate, while $N_{s}$ denotes the number of excited states with the same quantum numbers as the pion.

To investigate the asymptotic behavior of $C_{\text{2pt}}$ at large $t_{\rm sep}$, we define the effective mass as

$$m_{\rm eff}(t_{\rm sep})=\ln\left(\frac{C_{\text{2pt}}\left(t_{\rm sep}\right)}{C_{\text{2pt}}\left(t_{\rm sep}+a\right)}\right)~{}.$$ (23)

The effective masses for different momenta are shown in the upper panel of Fig. 1. As $t_{\rm sep}$ increases, the effective masses approach plateaus, which align with the dashed lines computed from the relativistic dispersion relation $E_{0}=\sqrt{m_{0}^{2}+(P^{z})^{2}}$. This agreement indicates that the ground state is effectively isolated from the excited-state tower for $t_{\rm sep}\gtrsim 8a$.

In practical calculations, the sum over $N_{s}$ must be truncated, as higher excited-state contributions decay rapidly with increasing $t_{\rm sep}$. In this work, we perform a two-state fit by setting $N_{s}^{\rm max}=2$, which allows us to efficiently extract the ground-state contribution while accounting for the leading excited-state contamination. The lower panel of Fig. 1 presents the extracted ground-state energies as a function of the hadron momentum. To test the validity of the relativistic dispersion relation, we fit the data points using the functional form

$$E=\sqrt{m_{0}^{2}+c_{1}P^{2}+c_{2}a^{2}P^{4}}~{},$$ (24)

where $E$ is the ground-state energy, $P$ is the hadron momentum, and $a$ is the lattice spacing. The coefficients $c_{1}$ and $c_{2}$ parameterize the possible discretization effects. As shown by the fit bands in Fig. 1, the extracted ground-state energies exhibit excellent agreement with the expected relativistic dispersion relation up to $P\approx 3.6$ GeV. This agreement shows that discretization effects remain small in the dispersion relation within the momentum range studied in this work.

IV.2 Bare quasi-TMD beam function matrix elements of pion

Figure 2: From left to right, the ratios of three-point to two-point functions $R_{\tilde{h}}$ for pion quasi-TMD with hadron momentum $P^{z}=1.83$, 2.43 and 3.04 GeV are shown as functions of $t_{\rm sep}$ and $\tau$. The upper and lower panels are for the cases with $(b_{\perp},z)=(1,3)a$ and $(b_{\perp},z)=(3,3)a$, respectively. The colored bands are two-state fit results while the gray band is the estimated ground-state matrix element.

Figure 3: The bare matrix elements of quasi-TMD in the coordinate space are plotted as the function of $b_{\perp}$ and $\lambda=zP^{z}$. These matrix elements are extracted from the two-state chained fit of $C_{\rm 2pt}$ and $R_{\tilde{h}}$. From left to right, the three panels correspond to hadron momenta of $P^{z}=1.83$, 2.43 and 3.04 GeV, respectively. It is observed that for all three hadron momenta, the quasi-TMD decays as the transverse separation $b_{\perp}$ increases and asymptotically approaches zero in the large $\lambda$ regime.

The quasi-TMD beam function is extracted from the three-point correlator

$$C_{\tilde{h}}\left(t_{\mathrm{sep}},\tau\right)=\left\langle\pi_{s}(\vec{x}_{0},t_{\rm sep})\bar{q}(z,b_{\perp},\tau)\Gamma q(\vec{0},\tau)\pi_{s}^{\dagger}(\vec{P},0)\right\rangle,$$ (25)

which can be expressed in terms of a spectral decomposition as

$$\begin{aligned} C_{\tilde{h}}\left(t_{\mathrm{sep}},\tau\right)=\sum^{N_{s}-1}_{n,m}\frac{z_{n}^{\dagger}O_{nm}z_{m}}{2\sqrt{E_{n}E_{m}}}e^{-E_{n}\left(t_{\text{sep}}-\tau\right)}e^{-E_{m}\tau}~{},\end{aligned}$$ (26)

where $O_{nm}=\langle n|\bar{q}(z,b_{\perp})\Gamma q(0)|m\rangle$ represents the matrix elements of the quasi-TMD operator. To extract the ground-state matrix element $O_{00}$, we employ a two-state chained fit using the two-point correlator $C_{\rm 2pt}$ and the ratio

$$R_{\tilde{h}}\left(t_{\mathrm{sep}},\tau\right)=\frac{C_{\tilde{h}}\left(t_{\mathrm{sep}},\tau\right)}{C_{\rm 2pt}\left(t_{\mathrm{sep}}\right)}~{}.$$ (27)

In the limit $t_{\rm sep}\to\infty$, the ratio $R_{\tilde{h}}$ asymptotically converges to the ground-state matrix element $O_{00}$. The chained fitting procedure— first fits the two-point correlator $C_{\rm 2pt}$, then uses the posterior distributions of $E_{0}$, $E_{1}$, $z_{0}$, and $z_{1}$ as priors in the subsequent fit of $R_{\tilde{h}}$.

For illustration, six examples of the chained fit applied to the quasi-TMD matrix elements are shown in Fig. 2. The error bars represent the data points of the ratio $R_{\tilde{h}}$, the colored bands depict the results of the two-state fit, and the gray band indicates the extracted ground-state matrix element. Additional details on the ground-state fit can be found in App. D.

The extracted bare matrix elements of the quasi-TMD function in the coordinate space are presented in Fig. 3. The transverse separation $b_{\perp}$ is plotted up to $1.2$ fm for three different hadron momenta: $P^{z}=1.83$, $2.43$, and $3.04$ GeV. In all cases, the quasi-TMD function decreases as $b_{\perp}$ increases and asymptotically approaches zero in the large-$z$ regime, consistent with expected physical behavior.

IV.3 Renormalization and extrapolation

Figure 4: The extrapolation of the renormalized quasi-TMD in coordinate space for various hadron momenta is shown in three panels. From left to right, they correspond to hadron momenta $P^{z}=1.83$, 2.43 and 3.04 GeV, respectively. The regions between the two red dashed lines indicate the extrapolation range using Eq. (31). For different hadron momenta, the same starting point of $z=0.78$ fm is chosen for the extrapolation.

As discussed in Ref. [77], the absence of Wilson lines in the CG correlator eliminates linear divergences, allowing the renormalized operator to be defined as

$$\bar{q}_{0}(z,b_{\perp})\Gamma q_{0}(0)=Z_{q}(a)\left[\bar{q}(z,b_{\perp})\Gamma q(0)\right]~{},$$ (28)

where $Z_{q}$ is the CG quark wave function renormalization factor. This renormalization is independent of both the external hadron states and the spatial separation of the quark bilinear operator. Consequently, an appropriate ratio can be constructed to cancel out the renormalization factor.

Following the approach in Ref. [66], we renormalize the quasi-TMD matrix elements using the ratio

$$\tilde{h}_{\gamma^{t}\gamma^{5}}(z,b_{\perp},P^{z};\mu)=\frac{\tilde{h}^{0}_{\gamma^{t}\gamma^{5}}(z,b_{\perp},P^{z};a)}{\tilde{\varphi}^{0}_{\gamma^{t}\gamma^{5}}(z=0,b_{\perp},P^{z}=0;a)}~{},$$ (29)

and similarly for the quasi-TMDWF,

$$\tilde{\varphi}_{\gamma^{z}\gamma^{5}}(z,b_{\perp},P^{z};\mu)=\frac{\tilde{\varphi}^{0}_{\gamma^{z}\gamma^{5}}(z,b_{\perp},P^{z};a)}{\tilde{\varphi}^{0}_{\gamma^{t}\gamma^{5}}(z=0,b_{\perp},P^{z}=0;a)}~{}.$$ (30)

Here, $\tilde{\varphi}^{0}$ denotes the bare quasi-TMDWF matrix elements defined in Eq. (11) and extracted in App. C.

The renormalized matrix elements in the coordinate space are presented in Fig. 4. As shown, the quasi-TMD matrix elements decay rapidly as a function of $\lambda=zP^{z}$, reaching approximately zero for $\lambda\gtrsim 5$. However, at large distances, while the values remain statistically consistent with zero, the statistical uncertainties persist at a constant level, which will lead to non-physical fluctuations in the direct Fourier transform. Due to the finite correlation length of spatial correlators in QCD [91], the quasi-TMD matrix elements in coordinate space are expected to exhibit exponential decay when the coordinate separation $z$ is large. Moreover, as demonstrated in Ref. [91], the extracted quasi-distributions in momentum space within the moderate $x$ region remain largely insensitive to the choice of extrapolation strategy. Therefore, we apply a non-fit extrapolation method to smooth the uncertainties of the renormalized quasi-TMD matrix elements at long distances. The extrapolation is performed using the following form:

$$\tilde{h}^{\rm ext}=w\cdot\tilde{h}+(1-w)\cdot 0~{},$$ (31)

where $w$ is a weight function that transitions linearly from 1 to 0 within the range $\lambda\in[z_{\rm ext}P^{z},z_{\rm max}P^{z}]$, the $z_{\rm ext}=0.78$ fm is the starting point of extrapolation and $z_{\rm max}=1.2$ fm is the largest longitudinal separation of quasi-TMD. It is expected that $z>z_{\rm ext}$ is large enough to see the exponential decay behavior of quasi-TMD. In addition, the comparison of TMDPDF using different $z_{\rm ext}$ can be found in App. E. The extrapolation range is indicated by the two red dashed lines in Fig. 4. After applying this extrapolation, the uncertainty bands smoothly converge to zero, mitigating non-physical fluctuations.

V Pion valence quark TMDPDF

As shown in the factorization formula Eq. (9), the computation of the unpolarized pion valence-quark TMDPDF relies on three key inputs: the Collins-Soper (CS) kernel, the intrinsic soft function, and the quasi-TMD beam function matrix elements discussed in Sec. IV. In this section, we present the numerical results for each of these components and, ultimately, the extracted unpolarized valence TMDWF and TMDPDF of the pion.

Figure 5: The ratio $\gamma^{\overline{\rm MS}}(b_{\perp},\mu,P_{1},P_{2})$ of different momentum pairs, for $b_{\perp}=2a$ (upper panel) and $b_{\perp}=8a$ (lower panel), at $\mu=2$ GeV.

V.1 The Collins-Soper kernel

The CS kernel can be extracted from the ratio of quasi-TMDWFs with different momenta, as described in Eq. (14). The bare matrix elements of the quasi-TMDWFs are extracted in App. C, which follows the same strategy as in Ref. [66]. For the matching procedure, we use the fixed-order one-loop results for the CG matching coefficient $C_{\text{TMD}}$ and the corresponding hard kernel $H_{\phi}$, as provided in Ref. [78]. Furthermore, we account for large logarithmic terms in the hard kernel by applying renormalization group evolution to improve the accuracy of the perturbative matching up to NLL; details on this resummation can be found in App. A.

Fig. 5 shows the CS kernel $\gamma^{\overline{\rm MS}}(b_{\perp},P_{1},P_{2};\mu)$ extracted from different momentum combinations $n_{1}^{z}/n_{2}^{z}$. The results are consistent across different momentum ratios within the uncertainty bands. However, a subtle $x$-dependence and variations between different momentum ratios suggest the presence of power corrections, since the pion mass $m_{\pi}=670$ MeV is quite heavy, contributing to systematic uncertainties. To estimate these uncertainties, we select two sets of closely spaced momentum values: $n^{z}_{1}/n^{z}_{2}=8/9$ and $n^{z}_{1}/n^{z}_{2}=9/10$. Within each Jackknife sample, we collect two momentum pairs and include 80 data points over the interval $x\in[0.34,0.66]$. The mean value $\langle\gamma\rangle_{i}$ and the standard deviation $\sigma_{i}$ are then calculated for each Jackknife sample $i$. The final mean value, along with the corresponding statistical and systematic uncertainties, is determined as follows:

$$\begin{aligned} \text{Mean}&=\text{Avg}[\langle\gamma\rangle_{i}]\\ \text{Stat}&=\text{Std}[\langle\gamma\rangle_{i}]\cdot\sqrt{N_{s}-1}\\ \text{Sys}&=\text{Avg}[\sigma_{i}]~{},\end{aligned}$$ (32)

where “Avg” means taking the average and “Std” means taking the standard deviation, the factor $\sqrt{N_{s}-1}$ is caused by the Jackknife resampling.

Figure 6: The CS kernel at $\mu=2$ GeV (red points with error bars) The 3-loop perturbative results [92, 93] for the CS kernel are denoted as $\rm N^{3}LO$ and $\rm N^{3}LL$ in the figure. The CS kernels from recent phenomenological parameterizations of experimental data are shown from MAP22 [19], IFY23 [20], ART23 [23], MAP24FI [25] and ART25 [26]. In addition, some recent lattice calculations are presented from LPC23 [63], ASWZ24 [65] and DWF24 [66].

After incorporating both statistical and systematic uncertainties, Fig. 6 presents our final results for the CS kernel at $\mu=2$ GeV, alongside comparisons with previous theoretical and phenomenological studies. In the small $b_{\perp}$ region, our results align well with the three-loop perturbative calculations from Refs. [92, 93], labeled $\rm N^{3}LO$ and $\rm N^{3}LL$. Moreover, our calculation remains reliable in the large $b_{\perp}$ region, where perturbative methods break down.

To further contextualize our findings, we compare them with CS kernels extracted from recent global fits of experimental data, including MAP22 [19], IFY23 [20], ART23 [23], MAP24FI [25] and ART25 [26]. Additionally, we include recent lattice QCD results from LPC23 [63], ASWZ24 [65], and DWF24 [66].

A particularly important observation is the strong agreement between our calculation and the DWF24 result from chirally symmetric domain-wall fermion configurations, both of which employ the CG framework despite using different lattice actions. Moreover, it is observable that the two most recent global analysis results (MAP24FI and ART25) exhibit a deviation from their prior results in different directions. Nevertheless, both remain consistent with this work due to the large uncertainty in our CS kernel, primarily stemming from power corrections at such a heavy pion mass, which resulted in non-flat curves in Fig. 5. In our future study, the precision of our CS kernel could be refined by adopting a smaller valence pion mass.

V.2 Intrinsic soft function

Figure 7: Extrapolations of reduced quasi-TMDWF using Eq. (33). From left to right, the three panels correspond to hadron momenta $P^{z}=3.44$, 3.87 and 4.30 GeV, respectively. All of three cases are evolved to $P_{\rm FF}=3.01$ GeV using the rapidity evolution of quasi-TMDWF. The shaded red bands indicate the regions constitute the input for the extrapolation fit.

Figure 8: Intrinsic soft function in CG at $\mu=2$ GeV in the ratio scheme, for different hadron momentum pairs of the quasi-TMDWF, $P_{\mathrm{WF}}$, and form factor, $P_{\mathrm{FF}}$. The final results are shown as pink points with error bars (see text for details). The corresponding one-loop fixed-order perturbative results are the red dashed line, and the one-loop RG-resummed (RGR) perturbative results are the blue dashed line.

The intrinsic soft function can be extracted from the analysis of form factors with large momentum transfer, aided by the quasi-TMDWF, as discussed in Sec. II.3. In this work, we use the form factors computed in Ref. [63], which include two datasets corresponding to the hadron momenta of $P_{\rm FF}=2.58$ GeV and $P_{\rm FF}=3.01$ GeV.

To extract the intrinsic soft function using Eq. (18), we integrate the quasi-TMDWF over $x_{1}$ and $x_{2}$. However, near the endpoint regions ($x\approx 0$ and $x\approx 1$), quasi-distributions suffer from significant power corrections. To mitigate this issue, we extrapolate the reduced quasi-TMDWF $\tilde{\Phi}$, using the functional form

$$\tilde{\Phi}(x)=cx^{n}(1-x)^{n}~{},$$ (33)

where $c$ and $n$ are fitting parameters. A comparison of the reduced quasi-TMDWF before and after extrapolation is shown in Fig. 7, where the three panels correspond to hadron momenta $P^{z}=3.44$, 3.87 and 4.30 GeV, respectively. All three cases are evolved to $P_{\rm FF}=3.01$ GeV using the rapidity evolution of quasi-TMDWF. The red shaded bands indicate the regions that constitute the input for the extrapolation fit, with the selected fit range of $x\in[0.3,0.45]$ and $x\in[0.55,0.7]$, where the results of three different momenta show consistency. The extrapolation form derived from the fitting process is applied to the endpoint regions with $x<0.3$ or $x>0.7$. Fig. 7 shows that the extrapolated results outside the fit range are in alignment with the observed trends. Furthermore, since the Sudakov kernel diverges near the endpoints ($x\approx 0$) after the RG resummation, the integral in Eq. (18) is restricted to the range $x_{1},x_{2}\in[0.05,0.95]$. This cutoff has a small impact on the results of the intrinsic soft function, since the integrand in the denominator of Eq. (18) converges near the endpoint regions.

To assess systematic effects from power corrections, we performed the integration with six different momentum pairs, formed by combining two momenta from the form factor dataset with three momenta from the quasi-TMDWF. In order to bridge the gap between $P_{\rm FF}$ and the hadron momenta of the quasi-TMDWF $P_{\rm WF}$, the quasi-TMDWF is evolved to $P_{\rm FF}$ by solving the rapidity evolution equation. The extracted intrinsic soft function, according to Eq. (18), is shown in Fig. 8, where comparisons across different momentum pairs indicate that power corrections have only a minor impact. Note that the superscript “ratio” of $S_{I}^{\rm ratio}$ indicates the ratio scheme in the renormalization of the quasi-TMDWF.

For the final results of $S_{I}^{\rm ratio}$, the central values and statistical uncertainties are determined from a correlated average over the six momentum pairs formed by combining the form factor and quasi-TMDWF datasets. The systematic uncertainties are determined by the spread of the mean values across these six momentum pairs, quantified as the absolute difference between the maximum and minimum values. The extracted intrinsic soft function at $\mu=2$ GeV, shown as pink points with error bars in Fig. 8, is compared to fixed-order and RG-resummed (RGR) perturbative predictions at leading-logarithmic (LL) accuracy. Further details on RG resummation can be found in App. A. In particular, at small $b_{T}$, our lattice results exhibit reasonable agreement with perturbative predictions, highlighting the robustness of our approach. In the large $b_{T}$ region, where perturbation theory becomes unreliable, our lattice results are expected to provide a more accurate description of the intrinsic soft function, offering valuable non-perturbative insights.

V.3 Pion TMDWF in $x$ space

Figure 9: The light-cone pion TMDWF at $y_{n}=0$, $P_{\mathrm{fix}}=4.30$ GeV and $\mu=2$ GeV of different hadron momenta, $P_{z}$, as the functions of momentum fraction, $x$, and for two transverse separations $b_{\perp}=4a$ (upper panel) and $b_{\perp}=8a$ (lower panel). The shaded gray bands ($x<0.11$ and $x>0.89$) indicate the endpoint regions where the estimated combined systematics are larger than $30\%$. The detailed estimation of the systematics is explained in App. F.

By integrating the CG kernel, intrinsic soft function, and the renormalized quasi-TMDWF, the light-cone TMDWF can be extracted by employing the factorization formula found in Eq. (12). Choosing scales $y_{n}=0$, $\zeta=(2xP_{\mathrm{fix}})^{2}$, and $\bar{\zeta}=(2\bar{x}P_{\mathrm{fix}})^{2}$, the factorization formula can be rewritten in a more explicit form,

$$\begin{aligned} &\sqrt{S_{I}\left(b_{\perp};\mu\right)}\cdot\tilde{\phi}_{\Gamma}\left(x,b_{\perp},P^{z};\mu\right)=\phi\left(x,b_{\perp};\mu,\zeta,\bar{\zeta}\right)\\ &\quad\quad\times H_{\phi}\left(x,\bar{x},P^{z};\mu\right)\exp\left[\ln\left(\frac{P^{z}}{P_{\mathrm{fix}}}\right)\gamma^{\overline{\mathrm{MS}}}\left(b_{\perp};\mu\right)\right]\\ &\quad\quad+\mathcal{O}\left(\frac{\Lambda_{\mathrm{QCD}}}{xP^{z}},\frac{1}{b_{\perp}\left(xP^{z}\right)},\frac{\Lambda_{\mathrm{QCD}}}{\bar{x}P^{z}},\frac{1}{b_{\perp}\left(\bar{x}P^{z}\right)}\right)~{},\end{aligned}$$ (34)

where the CS scale is evolved from $\zeta_{0}=(2xP^{z})^{2}$ and $\bar{\zeta}_{0}$ to $\zeta=(2xP_{\mathrm{fix}})^{2}$ and $\bar{\zeta}$ using the CS kernel extracted from quasi-TMDWFs.

The final results for the TMDWF at $y_{n}=0$, $P_{\mathrm{fix}}=4.30$ GeV, and $\mu=2$ GeV are shown in Fig. 9 as a function of the momentum fraction $x$ for three different hadron momenta. Two representative transverse separations, $b_{\perp}=4a$ (upper panel) and $b_{\perp}=8a$ (lower panel), are selected for illustration. As can be seen, the variation between different momenta remains mild in the moderate $x$ region, demonstrating the validity of power expansion in large $P^{z}$ within the quasi-TMD framework, where power corrections are small. The shaded gray bands ($x<0.11$ and $x>0.89$) indicate the endpoint regions where the estimated combined systematics are greater than $30\%$. The combined systematics include the variation of scales and power corrections, which are estimated by the variation of the initial scale $\mu_{0}$ of the RGR procedure by a factor of $\sqrt{2}$ and the spread in the central values of the TMDWFs with different momenta, respectively. More detailed discussion can be found in App. F.

V.4 Pion TMDPDF in $b_{\perp}$ and $k_{\perp}$ space

Figure 10: The unpolarized light-cone pion TMDPDF at $\zeta=\mu^{2}=4~{}\mathrm{GeV}^{2}$ of different hadron momenta are shown as the function of momentum fraction, $x$, and for transverse separations, $b_{\perp}=4a$ (upper panel) and $b_{\perp}=8a$ (lower panel). The shaded gray bands ($x<0.28$ and $x>0.81$) indicate the endpoint regions where the estimated combined systematics are larger than $30\%$. The detailed estimation of the systematics is explained in App. F.

Figure 11: Unpolarized light-cone TMDPDF of pion at $x=0.3$, $x=0.4$ and $x=0.5$ in $b_{\perp}$ space up to $b_{\perp}\gtrsim 1$ fm, and for differnt hadron momenta, $P_{z}$. The lattice results are illustrated alongside the phenomenological results from MAP22 [27] and JAM23 [28]. In addition, the collinear PDF from lattice calculation in CG (PDF Lat) [77] and global fit (PDF JAM21) [94] are plotted for comparison.

The combination of the CG kernel, intrinsic soft function, and renormalized quasi-TMD beam functions enables the extraction of the light-cone TMDPDF using the factorization formula in Eq. (9). The final results for the TMDPDF at $\zeta=\mu^{2}=4~{}\mathrm{GeV}^{2}$ are shown in Fig. 10 as a function of the momentum fraction $x$ for three different hadron momenta. Two representative transverse separations, $b_{\perp}=4a$ (upper panel) and $b_{\perp}=8a$ (lower panel), are selected for illustration. As can be seen, the variation between different momenta remains mild in the moderate $x$ region, demonstrating the validity of power expansion in large $P^{z}$ within the quasi-TMD framework, where power corrections are small. The shaded gray bands ($x<0.28$ and $x>0.81$) indicate the endpoint regions where the estimated combined systematics are greater than $30\%$, with details provided in App. F.

In Fig. 11, we examine the $b_{\perp}$ dependence of the TMDPDF at $x=0.3$, $x=0.4$, and $x=0.5$. The results for different hadron momenta exhibit good consistency, which is improved as the $x$ value approaches $x=0.5$, thereby confirming our expectation of the behavior of power corrections. For comparison, we include recent global fits of the pion TMDPDF from MAP22 [27] and JAM23 [28]. Although slight deviations in magnitude are observed, the general trend of the lattice results aligns well with the global fits in terms of $b_{\perp}$ dependence.

Since the collinear PDF and the $b_{\perp}\to 0$ limit of the unpolarized TMDPDF differ only by a perturbative expansion in $b_{\perp}$ [23], we also compare our results to collinear PDFs. Specifically, we include the collinear PDF extracted from lattice calculations in CG (PDF Lat) [77] and the global fit result from JAM21 (PDF JAM21) [94]. As shown in the plots, the collinear PDFs closely match the TMDPDF in the small $b_{\perp}$ region, supporting the expected theoretical behavior.

Figure 12: Unpolarized light-cone TMDPDF of pion at $x=0.3$, $x=0.4$ and $x=0.5$ in $b_{\perp}$ space after extrapolation to $b_{\perp}=3$ fm. The collinear PDF from lattice calculation in CG (PDF Lat) [77] is plotted at $b_{\perp}=$ 0 fm. The point at $b_{\perp}=0.06$ fm is obtained through cubic interpolation.

An important advantage of the CG method is the absence of linear divergence, which allows us to probe large $b_{\perp}$ regions with a reasonable SNR. At $b_{\perp}\gtrsim 1$ fm, the TMDPDF smoothly decays to values consistent with zero, facilitating a feasible Fourier transform into the $k_{\perp}$ space. In this study, we combined the collinear pion PDF result (PDF Lat) and the pion TMDPDF to interpolate the point at $b_{\perp}=a$, then applied a Gaussian form to extrapolate the large $b_{\perp}$ regime up to $3$ fm, the extrapolation form is

$$f(b_{\perp})=Ae^{-m\cdot b_{\perp}^{2}}~{},$$ (35)

where $A$ and $m$ are fit parameters. The extrapolated unpolarized light-cone TMDPDF of three momenta are plotted in Fig. 12.

Figure 13: Unpolarized light-cone TMDPDF of pion at $x=0.3$, $x=0.4$ and $x=0.5$ in $k_{\perp}$ space. The lattice results are illustrated alongside the phenomenological results from MAP22 [27] and JAM23 [28].

After extrapolation, we performed a discretized Fourier transform, with the results presented in Fig. 13. The general behavior of the TMDPDF extracted from the lattice remains consistent with global fits and decay as a function of $k_{\perp}$. The extrapolation form in Eq. 35 is inspired by the confinement in the transverse plane, while the model dependence of the results in the $k_{\perp}$ space will be investigated in future works with improved data precision.

VI Conclusion

In this work, we have presented the first lattice QCD calculation of the pion unpolarized valence-quark TMDPDF within the framework of LaMET. The calculation was performed on a 2+1 flavor QCD ensemble with a lattice spacing of $a=0.06$ fm and a pion mass of 300 MeV. By leveraging high statistics, off-axis momenta, a well-tuned boosted Gaussian smeared pion source, and the novel CG approach, we attained significant hadron momenta reaching 3 GeV. We computed the quasi-TMD beam function, extracted the CS kernel from quasi-TMDWF, and derived the intrinsic soft function. To enhance perturbative accuracy, we performed RGR of the Sudakov kernel at NLL order throughout our analysis. These advancements collectively enabled the determination of both the pion TMDWF and TMDPDF from the first principles.

Our results of the CS kernel and intrinsic soft function demonstrate consistency with perturbative calculations in the small-$b_{\perp}$ region, where perturbative methods remain valid. Additionally, the extracted CS kernel shows agreement with both previous lattice QCD studies and phenomenological fits of the experimental data.

By combining the renormalized quasi-TMDPDF, the CS kernel, and the intrinsic soft function, we obtained the pion valence TMDPDF across a range of $b_{\perp}$. Our results extend beyond $b_{\perp}\gtrsim 1$ fm, allowing for a feasible Fourier transform into momentum space. The extracted TMDPDFs in the $b_{\perp}$ and $k_{\perp}$ space exhibit qualitative agreement with phenomenological fits.

The outcome of this study highlights the efficacy of the CG quasi-TMD approach in probing the transverse momentum structure of hadrons. Our results provide a valuable first-principles complement to global fits, especially in the non-perturbative region where experimental data remain scarce. Note that the pion form factor used in this work was obtained from a different lattice setup; in future studies, we plan to compute it consistently within the same lattice ensemble to ensure full compatibility. Future improvements, such as the incorporation of finer lattice spacings, larger momenta, and improved control over systematic uncertainties, will further enhance the precision of lattice QCD determinations of TMDPDFs. Additionally, extending this methodology to other hadrons, including the nucleon, will provide deeper insights into the partonic structure of QCD bound states.

Appendix A Renormalization group resummation

A.1 RGR for TMD hard kernel

As mentioned in Sec. II, the quasi-TMD beam function and quasi-TMDWF can be matched to their light-cone counterparts using the TMD matching coefficient $C_{\mathrm{TMD}}$. In CG, the matching coefficient is calculated up to NLO fixed-order perturbation theory as  [78]

$$C_{\rm TMD}\left(xP^{z},\mu\right)=1+\frac{\alpha_{s}C_{F}}{4\pi}\left[-\frac{L_{p}^{2}}{2}-3L_{p}-12+\frac{7\pi^{2}}{12}\right]~{},$$ (36)

in which $L_{p}=\ln\frac{\mu^{2}}{(2xP^{z})^{2}}$ and $C_{F}=4/3$. The RG evolution (RGE) of the matching coefficient is given in the literature [98]

$$\frac{d}{d\ln\mu}\ln C_{\rm TMD}=\gamma_{H}=-\Gamma_{\rm cusp}L_{p}+\gamma_{C}~{}.$$ (37)

Combining with the fixed-order results in Eq. (36), we can solve to get the anomalous dimensions up to NLO

$$\Gamma_{\rm cusp}^{(1)}=2C_{F}\cdot\frac{\alpha_{s}(\mu)}{4\pi}~{},$$ (38)

and

$$\gamma_{C}^{(1)}=(-6C_{F})\cdot\frac{\alpha_{s}(\mu)}{4\pi}~{}.$$ (39)

Employing the defined QCD beta function $\beta$, one can incorporate the evolution of the coupling constant as described in $d\ln\mu=\frac{d\alpha_{s}}{\beta}$, and subsequently perform the integration as

$$\begin{aligned} &\ln C_{\rm TMD}(\mu)-\ln C_{\rm TMD}(\mu_{0})\\ &\quad=\int_{\alpha_{s}(\mu_{0})}^{\alpha_{s}(\mu)}\frac{d\alpha_{s}}{\beta(\alpha_{s})}\left[-2\Gamma_{\rm cusp}\int_{\alpha_{s}(\mu_{0})}^{\alpha_{s}}\frac{d\alpha^{\prime}_{s}}{\beta(\alpha^{\prime}_{s})}+\gamma_{C}(\alpha_{s})\right]~{},\end{aligned}$$ (40)

where we can choose the physical scale $\mu_{0}=2xP^{z}$. Due to the double logarithmic term $L^{2}_{p}$ in Eq. (36), the next-leading-log (NLL) result needs the tree-level $C_{\rm TMD}(\mu_{0})$, the one-loop $\gamma_{C}$ and the two-loop $\Gamma_{\rm cusp}$.

A.2 RGR for Sudakov kernel

The NLO fixed-order Sudakov kernel is given in Ref. [83, 99, 52]

$$C_{\rm Sud}(\mu)=1+\frac{\alpha_{s}C_{F}}{4\pi}\left[-\ln^{2}\frac{\mu^{2}}{Q^{2}}-3\ln\frac{\mu^{2}}{Q^{2}}-8+\frac{\pi^{2}}{6}\right]~{},$$ (41)

where $Q^{2}=-q^{2}=-(p_{1}-p_{2})^{2}$. The NLO RGE of Sudakov kernel can be derived as

$$\frac{d}{d\ln\mu}\ln C_{\rm Sud}(\mu)=\gamma_{1}(\mu)~{},$$ (42)

with the anomalous dimension

$$\gamma_{1}^{(1)}(\mu)=-\frac{\alpha_{s}(\mu)C_{F}}{4\pi}\left[4\ln\frac{\mu^{2}}{Q^{2}}+6\right]~{}.$$ (43)

Similarly, incorporating the evolution of coupling constant via $d\ln\mu=\frac{d\alpha_{s}}{\beta}$, one can do the integral as

$$\begin{aligned} &\ln C_{\rm Sud}(\mu)-\ln C_{\rm Sud}(\mu_{0})\\ &\quad=\int_{\alpha_{s}(\mu_{0})}^{\alpha_{s}(\mu)}\frac{d\alpha_{s}}{\beta(\alpha_{s})}\left[-\frac{\alpha_{s}C_{F}}{4\pi}\left(8\int_{\alpha_{s}(\mu_{0})}^{\alpha_{s}}\frac{d\alpha^{\prime}_{s}}{\beta(\alpha^{\prime}_{s})}+6\right)\right]~{},\end{aligned}$$ (44)

where we can choose the physical scale $\mu_{0}=Q$. Due to the double logarithmic term in Eq. (41), the NLL result needs the tree-level $C_{\rm Sud}(\mu_{0})$ and the one-loop $\gamma_{1}$.

A.3 RGR for intrinsic soft function

In $\overline{\rm MS}$ scheme, the NLO fixed-order intrinsic soft function in CG is [78]

$$S_{I}(b_{\perp},\mu)=1-\frac{\alpha_{s}C_{F}}{\pi}L_{b}~{},$$ (45)

where $L_{b}=\ln\frac{\mu^{2}b_{\perp}^{2}e^{2\gamma_{E}}}{4}$. The RGE can be derived as

$$\frac{d}{d\ln\mu}\ln S_{I}(b_{\perp},\mu)=-2\frac{\alpha_{s}(\mu)C_{F}}{\pi}~{}.$$ (46)

Similarly, incorporating the evolution of coupling constant via $d\ln\mu=\frac{d\alpha_{s}}{\beta}$, one can do the integral as

$$\begin{aligned} &\ln S_{I}(b_{\perp},\mu)-\ln S_{I}(b_{\perp},\mu_{0})=-2\frac{C_{F}}{\pi}\int_{\alpha_{s}(\mu_{0})}^{\alpha_{s}(\mu)}\frac{d\alpha_{s}}{\beta}\alpha_{s}~{}.\end{aligned}$$ (47)

The leading-log (LL) result is given with the combination of the tree-level $S_{I}(b_{\perp},\mu_{0})$ and the one-loop anomalous dimension.

Appendix B Scheme conversion

As discussed in Sec. II, using different renormalization schemes for quasi-TMDWF, one can get the intrinsic soft function in the corresponding renormalization schemes. This section derives the scheme conversion between the $\overline{\mathrm{MS}}$ scheme and the ratio scheme as defined in renormalization Eqs. (29) and (30).

Using the ratio scheme, the renormalized quasi-TMDWF is defined as

$$\tilde{\varphi}^{\text{ratio }}\left(z,b_{\perp},P^{z};\mu\right)=\frac{\tilde{\varphi}^{0}\left(z,b_{\perp},P^{z};a\right)}{\tilde{\varphi}^{0}\left(z=0,b_{\perp},P^{z}=0;a\right)}=\frac{\tilde{\varphi}^{\overline{\rm MS}}\left(z,b_{\perp},P_{z};\mu\right)}{\tilde{\varphi}^{\overline{\rm MS}}\left(z=0,b_{\perp},P^{z}=0;\mu\right)}~{},$$ (48)

$$\tilde{\phi}^{\text{ratio }}(x,b_{\perp},P^{z};\mu)=P^{z}\int\frac{dz}{2\pi}e^{iz(xP^{z})}\tilde{\varphi}^{\text{ratio }}(z,b_{\perp},P^{z};\mu)=\int\frac{d\lambda}{2\pi}e^{i\lambda x}\frac{\tilde{\varphi}^{\overline{\rm MS}}\left(z,b_{\perp},P_{z};\mu\right)}{\tilde{\varphi}^{\overline{\rm MS}}\left(z=0,b_{\perp},P^{z}=0;\mu\right)}~{},$$ (49)

in which the denominator can be approximately replaced by the short distance factorization coefficient $C_{0}$ in CG, then we got

$$\tilde{\phi}^{\text{ratio }}(x,b_{\perp},P^{z};\mu)=\frac{\tilde{\phi}^{\overline{\rm MS}}(x,b_{\perp},P^{z};\mu)}{C_{0}(b_{\perp};\mu)}~{}.$$ (50)

Combining with the Eq. (18), we have the scheme conversion of intrinsic soft function as

$$S_{I}^{\rm ratio}(b_{\perp};\mu)=S_{I}^{\overline{\rm MS}}(b_{\perp};\mu)\cdot C^{2}_{0}(b_{\perp};\mu)~{}.$$ (51)

The NLO fixed-order result of $C_{0}$ can be found in Ref. [77] as

$$C_{0}\left(b_{\perp};\mu\right)=1+\frac{\alpha_{s}C_{F}}{2\pi}\left(\frac{1}{2}-\frac{L_{b}}{2}\right)~{}.$$ (52)

One can do the RGR in the similar way as in App. A to get the LL result of $C_{0}$ as

$$C_{0}\left(b_{\perp};\mu\right)=\exp\left(\frac{C_{F}}{\beta_{0}}\ln\frac{\alpha_{s}(\mu)}{\alpha_{s}\left(2e^{-\gamma_{E}}/b_{\perp}\right)}\right)~{}.$$ (53)

Appendix C Bare matrix elements of quasi-TMDWF

Figure 14: The ground-state energies $E_{0}$ extracted from two-state fit of the two-point functions. The red line represents the exact dispersion relations with the valence pion mass $m_{0}=670$ MeV.

Figure 15: From left to right, the ratios of quasi-TMDWF to two-point functions $R_{\tilde{\varphi}}$ with hadron momentum $P^{z}=3.44$, 3.87 and 4.30 GeV are shown as functions of $t_{\rm sep}$. The upper and lower panels are for the cases with with $(b_{\perp},z)=(1,4)~{}a$ and $(b_{\perp},z)=(3,5)~{}a$, repectively. The colored bands are joint fit results.

Figure 16: Bare matrix elements of quasi-TMDWF in the coordinate space. The non-zero-momentum correlators extracted from two-state joint fits of $C_{\text{2pt}}$ and $C_{\tilde{\varphi}}$, while the zero-momentum correlators are extracted from one-state fit of $C_{\tilde{\varphi}}$. All bare matrix elements are normalized using the mean value of the local matrix element.

Figure 17: Quasi-TMDWF of pion in the momentum space. From up to down, the hadron momenta are $P^{z}$ = 3.44, 3.87 and 4.30 GeV. The small bumps in the lowest subfigure are caused by the small non-zero tails in the coordinate space, primarily manifest at large $b_{\perp}$.

The quasi-TMDWF is extracted from the non-local two-point correlator

$$C_{\tilde{\varphi}}\left(t_{\rm sep}\right)=\left\langle\bar{q}(z,b_{\perp},t_{\rm sep})\Gamma q(\vec{0},t_{\rm sep})\pi_{s}^{\dagger}(\vec{P},0)\right\rangle,$$ (54)

which can be expressed in terms of a spectral decomposition as

$$\begin{aligned} &C_{\tilde{\varphi}}\left(t_{\rm sep}\right)\\ &\quad=\sum_{n=0}^{N_{\mathrm{s}}-1}\frac{z_{n}}{2E_{n}}\langle\Omega|O^{\rm CG}_{\Gamma}|n\rangle\left(e^{-E_{n}t_{\rm sep}}+e^{-E_{n}\left(L_{t}-t_{\rm sep}\right)}\right)~{},\end{aligned}$$ (55)

where $\ket{\Omega}$ is the interactive vacuum state, $\ket{n}$ is the $n$-th energy eigenstate. When $\Gamma=\gamma^{z}\gamma^{5}$, the matrix element of the ground state is $\langle\Omega|O^{\rm CG}_{\gamma^{z}\gamma^{5}}|0\rangle=if_{\pi}P^{z}\tilde{\varphi}^{0}_{\gamma^{z}\gamma^{5}}$, where $f_{\pi}$ is the bare pion decay constant, and $\tilde{\varphi}^{0}_{\Gamma}$ is the bare matrix element of quasi-TMDWF in the coordinate space. When $\Gamma=\gamma^{t}\gamma^{5}$, the matrix element of the ground state is $\langle\Omega|O^{\rm CG}_{\gamma^{t}\gamma^{5}}|0\rangle=f_{\pi}E_{0}\tilde{\varphi}^{0}_{\gamma^{t}\gamma^{5}}$, where $E_{0}=\sqrt{m_{\pi}^{2}+(P^{z})^{2}}$.

In order to determine the energy spectrum and facilitate the extraction of the quasi-TMDWF, the local two-point functions of the pion are additionally computed with a valence pion mass of $m_{\pi}=670$ MeV for the corresponding hadron momenta. The ground-state energies from the two-state fit of the two-point functions are illustrated in Fig. 14, in comparison with the exact dispersion relations with $m_{0}=670$ MeV. Good consistency is found up to the hadron momentum $P^{z}=4.30$ GeV.

For correlators with nonzero momenta, we employ the two-state joint fit of $C_{\text{2pt}}$ and $C_{\tilde{\varphi}}$, while for zero-momentum correlators, the one-state fit is adopted to fit $C_{\tilde{\varphi}}$. Some of the two-state fits are selected as examples in Fig. 15. As a means of illustration, both data and fit results are plotted as the ratio of $C_{\tilde{\varphi}}$ and $C_{\text{2pt}}$, which is defined as

$$R_{\tilde{\varphi}}\left(t_{\rm sep}\right)=\frac{C_{\tilde{\varphi}}\left(t_{\rm sep}\right)}{C_{\text{2pt}}\left(t_{\rm sep}\right)}~{}.$$ (56)

This ratio will converge to $\langle\Omega|O^{\rm CG}_{\Gamma}|0\rangle/z_{0}$ in the $t_{\rm sep}\to\infty$ limit. More details on the ground-state fit can be found in App. D.

The bare matrix element of quasi-TMDWF in the coordinate space can be found in Fig. 16, different momenta are normalized using the mean value of the local matrix element at $(z,b_{\perp})=(0,0)$, so that some artifacts such as discretization effects can be canceled out. Due to the good convergence in the large-$\lambda$ region, the renormalized matrix elements $\varphi_{\Gamma}$ can be Fourier transformed directly to the momentum space.

Fig. 17 illustrates the quasi-TMDWF in momentum space, where it can be seen that the dependence on momentum is not significant, implying that the CS kernel is not a large quantity. The small bumps in the lowest subfigure are caused by the small non-zero tails in the coordinate space, primarily manifest at large $b_{\perp}$.

Appendix D Ground state fit

Figure 18: The evaluation of joint fits of quasi-TMDWF matrix elements. The left panel is the density distribution of $\chi^{2}/\rm{d.o.f.}$ and the right panel is the cumulative distribution function (CDF) of p-value.

Figure 19: The evaluation of chained fits on $R_{\tilde{h}}$ of quasi-TMD matrix elements. The left panel is the density distribution of $\chi^{2}/\rm{d.o.f.}$ and the right panel is the cumulative distribution function (CDF) of p-value.

We have performed a fully correlated Bayesian analysis of the two-point and three-point correlation functions to extract the bare matrix element $\tilde{h}^{0}_{\Gamma}(x,b_{\perp},P^{z};\mu)$ defined in Eq. (2) and $\tilde{\varphi}^{0}_{\Gamma}(x,b_{\perp},P^{z};\mu)$ defined in Eq. (10). The correlation across different data sets is taken into account by performing a Bayesian least-squares fit on each sample of Jackknife resampling. The parameter settings for the ground state fit are collected in Table 1.

Fit	Parts	$N_{s}$	$t_{\rm sep}$ range	$\tau$ range
$\tilde{\varphi}^{0}_{\gamma^{t}\gamma^{5}}$	$C_{\tilde{\varphi}}(t_{\rm sep})$	1	$t_{\rm sep}\in[12,15]$	/
$\tilde{\varphi}^{0}_{\gamma^{z}\gamma^{5}}$	$C_{\rm 2pt}(t_{\rm sep})$	2	$t_{\rm sep}\in[3,10]$	/
	$C_{\tilde{\varphi}}(t_{\rm sep})$	2	$t_{\rm sep}\in[5,9]$	/
$\tilde{h}^{0}_{\gamma^{t}\gamma^{5}}$	$C_{\rm 2pt}(t_{\rm sep})$	2	$t_{\rm sep}\in[3,12]$	/
	$R_{\tilde{h}}(t_{\rm sep},\tau)$	2	$t_{\rm sep}\in\{6,8,10,12\}$	$\tau\in[3,t_{\rm sep}-3]$

Table 1: Collection of ground state fit settings. $N_{s}$ = 1 means that only one ground state is included in the fit functions. The $t_{\rm sep}$ range of $C_{\tilde{\varphi}}$ varies in some sets of $(P^{z},b_{\perp})$ due to the different behaviors of excited states.

To evaluate the overall quality of the ground state fits, the density distribution of $\chi^{2}/\rm{d.o.f.}$ and the cumulative distribution function (CDF) of the p-value are plotted in Figs. 18 and 19, corresponding to quasi-TMDWF and qTMDPDF, respectively. The figures reveal that $\chi^{2}/\rm{d.o.f.}$ is predominantly distributed within the expected range, with p-values exceeding $0.05$ for the majority of fits, thereby suggesting a high level of overall quality. In addition, as examples, the fit results of the quasi-TMDWF and the quasi-TMD are plotted alongside the data points in Figs. 2 and 15, respectively.

Appendix E Stability of extrapolation

Figure 20: The unpolarized light-cone pion TMDPDF of different $z_{\rm ext}$ are ploted in comparison as the function of momentum fraction $x$. Two hadron momenta $P^{z}=2.43$ GeV (upper panel) and $P^{z}=3.04$ GeV (lower panel) are selected as examples. The variations among different $z_{\rm ext}$ are mild in the moderate region near $x=0.5$. The shaded gray bands ($x<0.28$ and $x>0.81$) indicates the endpoint regions where the estimated combined systematics are larger than $30\%$. The detailed estimation of the systematics is explained in App. F.

As mentioned in Sec. IV, the light-cone TMDPDF in momentum space within the moderate $x$ region is insensitive to extrapolation strategies. The unpolarized light-cone pion TMDPDF with different $z_{\rm ext}$ are plotted in Fig. 20 for comparison.

It shows that the variations of the TMDPDF among the different $z_{\rm ext}$ are mild in the moderate region near $x=0.5$, which shows the stability of the extrapolation.

Appendix F Estimation of systematics

Figure 21: The estimation of the relative systematics of TMDWF (left panel) and TMDPDF (right panel). The scale variation is estimated by varying the initial scale $\mu_{0}$ of the RGR procedure by a factor of $\sqrt{2}$. The momentum variation is quantified as the ratio of the spread in the central values of the light-cone distributions, evaluated at three different momenta, to their mean. The combined systematics are estimated by the root-sum-square of two variations. A threshold of $30\%$ on the combined systematics is used to define the reliable region, yielding $x\in[0.11,0.89]$ for the TMDWF and $x\in[0.28,0.81]$ for the TMDPDF. These regions are indicated by the red markers in the figure.

As mentioned in Sec. II, the results of TMDWF and TMDPDF are only reliable in the moderate $x$ region due to various systematics such as uncontrolled power corrections in the endpoint regions. In this section, we will give an approximate estimation of the systematics based on the results of light-cone distributions, with the objective of delineating the reliable range within the $x$-space.

The estimation of the systematics is separated into two parts, the systematics of varying the initial scale of RGR, and the variation of light-cone distributions calculated at different pion momenta. Taking TMDWF and TMDPDF at $b_{\perp}=6a$ as representatives, the relative systematics are plotted in Fig. 21. The systematics of scale variation is estimated by varying the initial scale $\mu_{0}$ of the RGR procedure by a factor of $\sqrt{2}$. Note that the systematics of scale variation is dependent on the hadron momentum. For estimation purposes, the maximum hadron momentum is selected as the representative. The momentum variation is quantified as the ratio of the spread in the central values of the light-cone distributions, evaluated at three different momenta, to their mean. These two sources of systematic uncertainty are regarded as independent, so they are combined using the root-sum-square method to provide an estimation of the combined systematics. A threshold of $30\%$ on the combined systematics is used to define the reliable region, yielding $x\in[0.11,0.89]$ for the TMDWF and $x\in[0.28,0.81]$ for the TMDPDF.