Nonperturbative Collins-Soper Kernel from Chiral Quarks with Physical Masses

We present a lattice QCD calculation of the rapidity anomalous dimension of quark transverse-momentum-dependent distributions, i.e., the Collins-Soper (CS) kernel, up to transverse separations of about 1 fm. This unitary lattice calculation is conducted, for the first time, employing the chiral-symmetry-preserving domain wall fermion discretization and physical values of light and strange quark masses. The CS kernel is extracted from the ratios of pion quasi-transverse-momentum-dependent wave functions (quasi-TMDWFs) at next-to-leading logarithmic perturbative accuracy. Also for the first time, we utilize the recently proposed Coulomb-gauge-fixed quasi-TMDWF correlator without a Wilson line. We observe significantly slower signal decay with increasing quark separations compared to the established gauge-invariant method with a staple-shaped Wilson line. This enables us to determine the CS kernel at large nonperturbative transverse separations and find its near-linear dependence on the latter. Our result is consistent with the recent lattice calculation using gauge-invariant quasi-TMDWFs, and agrees with various recent phenomenological parametrizations of experimental data.


Introduction
The parton transverse-momentum-dependent distributions (TMDs) are crucial for a three-dimensional understanding of parton motions within a hadron, offering a more comprehensive view than traditional onedimensional parton distribution functions (PDFs).It sheds light not only on the intrinsic motion of partons in transverse directions but also on the interplay between the transverse momentum of quarks and the spin of nucleons or quarks themselves.This comprehensive perspective is crucial for a deep understanding of the dynamic and complex nature of nucleons.The accurate characterization of TMDs is also critical for interpreting experimental data from high-energy collisions, particularly in relation to the transverse momentum distributions of electroweak and Email address: gaox@anl.gov(Xiang Gao) Higgs bosons [1,2].They are fundamental to precision measurements, such as determining the mass and width of the W boson [3,4].As high-energy physics experiments continue to advance, the measurement of TMDs will become increasingly important.The ongoing and future experiments at facilities such as the Large Hadron Collider [5,6] and Electron-Ion Collider [7,8,9,10] are expected to profoundly enrich our knowledge of TMDs.This will not only enhance our grasp of hadron structure and nucleon spin but also contribute significantly to the broader field of particle physics.
Central to the study and practical application of TMDs is the Collins-Soper (CS) kernel, which is responsible for the (rapidity) scale evolution of TMDs [11,12], enabling the consistent interpretation of experimental data across different energy scales.It is instrumental in connecting theoretical predictions with experimental observations.The TMDs and CS kernel can be extracted through the global analysis of experimental data including the semiinclusive deep inelastic scattering (SIDIS) and Drell-Yan processes [13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29].However, the nonperturbative nature of Quantum Chromodynamics (QCD) at low transverse momenta necessitates certain parametrizations and introduces model dependence.As a result, there is an increasing interest for these intrinsically nonperturbative quantities to be calculated directly from first-principles lattice QCD.
Although direct simulation of TMDs on the Euclidean lattice is impractical, as TMDs are defined on the lightcone, it has been demonstrated that they can be accessed through quasi-TMDs within the framework of Large-Momentum Effective Theory (LaMET) [30,31,32].The quasi-TMDs involve the matrix elements of equal-time gauge-invariant (GI) operators: with b = (b ⊥ , b z ) covering both longitudinal (b z ) and transverse (b ⊥ ) directions, linked by a staple-shaped Wilson line W ⊐ whose length is characterized by η.In the large momentum and η → ∞ limit, the quasi-TMDs can be related to the light-cone TMDs through the perturbative factorization [33,34,35,36,37,36,38,39,40,41,42,43,44,45,46].Building on this, significant advancements have been achieved over the past few years.The CS kernel has been extracted from either quasi-TMD parton distribution functions (TMD-PDFs) [47,48,49], quasi-TMD wave functions (TMD-WFs) [50,51,52,49,53,54,55], or the moments of the quasi-TMDs [56].Lattice QCD calculations of soft functions [50,51,53], along with the first results of nucleon TMD PDFs [57] and pion TMDWFs [58] also have been reported.Additionally, progress has also been made in the systematical control of these calculations, including improved matching up to two loops [59,60], addressing the operator mixing and working at physical quark masses [61,62,63,64,65,66,54,67,55].Despite notable progress, the lattice calculation of TMDs remains challenging.To suppress power corrections, a large momentum P z is required, which incurs a significant computational cost.In addition, the signalto-noise ratio of the quasi-TMD matrix elements is adversely affected by exponential decay as the total length of the space-like Wilson line increases.This decay makes it particularly difficult to investigate quasi-TMDs at large b ⊥ , where the results are desired to complement the phenomenological analysis.What's more, the linear divergence and pinch-pole singularity in the Wilson lines also complicate the renormalization procedure, although they could be cancelled by the Wilson loop [68] or held fixed by keeping a constant length of the Wilson line for given b ⊥ [54].Besides, the operator mixings of Wilson-line operators [62,63,64,65,66] also need to be subtracted systematically [54].Recently, a novel approach has been proposed for computing parton physics in the Coulomb gauge (CG) [69], notably without the use of Wilson lines.Thereby, the complexity induced by the Wilson line can be avoided.It has been demonstrated that, in the large momentum limit, the CG quasi-PDF falls into the same universality class as the GI case under the LaMET framework.Further progress has been made in the realm of quasi-TMDs [70], involving the matrix elements of equaltime operators, with the CG condition ∇•A = 0 but without a Wilson line.
The factorization of quasi-TMDs in the CG has been derived from the soft collinear effective theory (SCET) [70] and verified at one loop in perturbation theory [70,71].
In this study, we have, for the first time, computed the quasi-TMDWFs of the pion in the CG and extracted the CS kernel from these measurements.Without Wilson lines, the CG correlators are multiplicatively renormalizable and free from the linear divergence [69] and pinch singularity, as well as the operator mixings originating from the Wilson line geometry.Through our calculation, we show that, the CG approach leads to consistent CS kernel with the conventional GI approach.Moreover, the CG approach can significantly reduce the signal-to-noise ratio and extend the prediction power of lattice computation in the nonperturbative regime of interest to TMD physics.

Theoretical framework
The pion quasi-TMDWF in the CG is defined as the Fourier transform of the matrix elements: where the pion is boosted with momentum P = (0, 0, P z ).By selecting Γ = γ t γ 5 or γ z γ 5 , the quasi-TMDWF φΓ (x, b ⊥ , P z , µ) can be related to the light-cone TMDWF ϕ(x, b ⊥ , ζ, µ) (under the principle-value prescription of the light-cone Wilson lines [72]) in the large P z limit through perturbative factorization, which can be expressed as [37,42,70], ) is a hard matching kernel that has been computed from one-loop perturbation theory [70,71].S r (b ⊥ , µ) represents the reduced soft functions, extractable from the form factors of fast-moving color-charged states [37,42,70].Consequently, the xdependent light-cone TMDWF can be derived, subject to power corrections that are suppressed by large P z and b ⊥ .Alternatively, the CS kernel can be extracted through the ratios of the quasi-TMDWFs with different momenta P 1 and P 2 [33,35,50], with perturbative corrections δγ MS inferred from H(x, x, P z , µ) and power corrections (p.c.).

Lattice setup
The bare matrix elements of the pion quasi-TMDWF can be extracted from the two-point correlation functions in the lattice simulations.For CG quasi-TMDWFs, we compute, O CG Γ (b, P, t s ) Here y 0 is the source position, and t s is the time separation.We chose Γ = γ t γ 5 , as it should be free from the operator mixings caused by chiral symmetry breaking [62,64,65] under our lattice setup.
To improve the signal-to-noise ratio and increase the overlap with the pion ground state, we used extended pion source after boosted Gaussian smearing [73], with the s denoting the smeared fields.We also compute the pion-pion two-point functions, with smeared source and sink to extract the energy spectrum created by π † as well as the overlap amplitudes.
For the lattice simulation, we utilized a 2+1-flavor Domain-wall gauge ensemble generated by RBC and UKQCD Collaborations of size N 3 s ×N t ×N 5 = 64 3 ×128× 12, denoted by 64I [74].The quark masses are at the physical point and the lattice spacing is a −1 = 2.3549 (49) GeV (a = 0.0836 fm).For the boosted Gaussian smearing, the Gaussian radius was chosen to be r G = 0.58 fm, and we chose the quark boost parameter j z to be 0 and 6 [75,76] which are optimal to hadron momentum P z = 2πn z /(N s a) with n z = 0 and 8, so that the largest momentum in our calculation is P z = 1.85 GeV.Since only two-point functions are involved in this calculation, measurements at other momenta (n z ∈ [0, 8]) were also computed through contractions using the same profiled quark propagator.
To increase the statistics, we used 64 configurations coupled with All Mode Averaging (AMA) technique [77].We computed 2 exact and 128 sloppy solutions for the quasi-TMDWFs with momenta n z ∈ [4,8], while 1 exact and 32 sloppy solutions for the cases with n z = [0, 3].The quark propagators are evaluated from CG-fixed configurations using deflation based solver with 2000 eigen vectors.
After fixing the CG, the GI quasi-TMDWF defined from Eq. ( 1) shares the same quark propagators as the CG case but needs an additional staple-shaped Wilson line to maintain the gauge invariance.Therefore, we also computed the GI quasi-TMDWF correlators C GI πO (t s ; b ⊥ , b z , P z , η) during the contraction.We chose η = 12a in this case using the same setup of the stapleshaped Wilson line as Ref. [54].We employed Wilson flow [78], with a flow time t F = 1.0 (roughly smears the gauge fields over the radius √ 8a 2 ), to suppress the ultraviolet (UV) fluctuations and enhance the signal-to-noise ratio.

Quasi-TMDWF
The pion-pion and quasi-TMDWF correlators have the following spectral decompositions, where E n (P z ) is the energy level, and Z n = ⟨n|π † (P z )|Ω⟩ is the overlap amplitude created by the pion interpolator (real and positive [54]).|Ω⟩ represents the vacuum state, while |n⟩ = |0⟩, |1⟩, ... represents the ground state as well as the excited states.To take the advantage of high correlations between the pion-pion and quasi-TMDWF two-point functions, we construct their ratio as, In Fig. 1, the ratios of our largest momentum n z = 8 at b ⊥ = 10a are shown as an example.In the t s → ∞ limit, this ratio will reduce to ⟨Ω|O γ t γ 5 |0⟩/Z 0 = E 0 φB /Z 0 and gives the bare quasi-TMDWF matrix elements φB (b ⊥ , b z , P z , a).In practice, with finite t s we truncate Eq. ( 10) and Eq. ( 11) up to N st = 2 and extract the bare matrix elements through the two-state fit.The fit results are shown as the bands in Fig. 1 which can nicely describe the data point.
The renormalization of CG quasi-TMD operator is straightforward.It involves only the CG quark wave function renormalization, which is an overall multiplicative constant and does not depend on the spatial separations b ⊥ and b z [69].
In contrast, the GI quasi-TMD operator defined in Eq. ( 1), though also subject to multiplicative renormalization, requires the removal of pinch pole singularities, cusp divergences, and linear divergences associated with the Wilson line [61,79,80,68].This renormalization is proportional to the total link length.In our implementation, the total length of the staple-shape Wilson line is 2η + b ⊥ , independent of b z .
Since the renormalization process solely involves the UV properties of operators and is independent of the external hadron states, we use the renormalization group invariant ratios [54] Φ without affecting the x-and P z −dependences of the quasi-TMDWF after Fourier transform.The above ratio also may reduce some correlated uncertainties and eliminate some power corrections.Thus, we also adopt the same procedure for the CG matrix elements, whose renormalization is b ⊥ and η independent.In Fig. 2, we show the renormalized matrix elements for our largest momentum, P z = 1.85 GeV, and for b ⊥ = 2a, 6a and 8a, as a function of b z , for both CG (filled squared symbols) and GI (open circled symbols) cases.) The left panel and right panel show the real and imaginary parts, respectively.It is evident that reasonable signal remains for the CG case even when b ⊥ become large.
In contrast, the signal-to-noise ratio of the GI matrix elements rapidly deteriorates as b ⊥ increases, primarily due to the long Wilson line and its UV fluctuations.In addition, it is shown that the imaginary parts of the CG case are consistently zero, while they have non-zero values for GI case.This is expected as the imaginary part depends on the longitudinal orientation of the Wilson line in the GI case [72], whereas the CG condition does not favor any direction [70].
One can also observe that the matrix elements decrease as a function of b z , diminishing to zero within the errors when b z ≳ 1 fm.This behavior facilitates the numerical Fourier transform to x-space with a simple truncation at the maximum value of b z , which is expressed as, where we apply a first-order spline interpolation to smooth the data points.Since the CG quasi-TMDWF correlator is real and symmetric in b z , the distribution must be real in the x-space.
In the upper panel of Fig. 3, we show selected results of the CG quasi-TMDWFs with momentum n z ∈ [4,8] and b ⊥ = 2a, 8a, 10a.Encouragingly, reasonable signal persists even when b ⊥ is as large as 10a.However, the signalto-noise ratio decreases as the momentum increases.In addition, it is evident that the quasi-TMDWFs, though ap-pearing to be non-zero outside the physical region, have a trend to shrink into x ∈ [0, 1] as the momentum increases.This observation is consistent with the power expansion of the LaMET, suggesting the quasi-TMDWF is approaching the light-cone TMDWF in the large momentum limit.

The Collins-Soper kernel
According to Eq. ( 5), we define the following estimator of the CS kernel utilizing the quasi-TMDWFs at finite momenta, In this work, we applied the perturbative corrections δγ MS derived from the next-to-leading logarithm (NLL) matching kernels for the CG case [38,43,54,70] as only oneloop non-cusp anomalous dimension is available.The MS scale has been set to be µ = 2 GeV.If the power corrections and higher-order corrections are small, γMS should be independent of P z and x.
In the lower panels of Fig. 3, we show the CS kernel estimators for the CG case with various combination of momenta, n 1 and n 2 , as a function of x.The x-independent plateaus can be found in the moderate x region within the  errors, which is robust even at the largest b ⊥ .This indicates the effectiveness of the factorization formula in (4).In the end-point regions of both small and large x, the results appear to diverge, signaling a breakdown of the factorization in these areas.However, the length of plateaus extend as the momentum increases, which is consistent with the power corrections suggested in Eq. ( 4).
As for the momentum dependence, it is absent for the case of large b ⊥ , which indicates well suppressed power corrections by 1/(P z b ⊥ ), despite their slightly larger errors.However, results at small b ⊥ (e.g., for 2a) and with small momentum (e.g., for n 1 = 4) deviated from the ones derived from larger momenta.This momentum dependence is reduced when n 1 and n 2 gets close, and disappear when n 1 and n 2 are close enough (e.g., for n 1 /n 2 = 4/5).This observation suggests that the power corrections and higher-order perturbative corrections not well suppressed in the cases of small b ⊥ and large differences in P .
To estimate the CS kernel, we averaged over the estimator γMS (x, b ⊥ , µ, P 1 , P 2 ) within x ∈ [x 0 , 1 − x 0 ] across various n 1 and n 2 .Only the cases of n 2 − n 1 = 1 considered.The value of x 0 is determined by requiring 2x 0 P z b ⊥ > 1 and 2x 0 P z > 0.7 GeV, suggested by the power correction.As a result, b ⊥ = a is always excluded in this work.The averages over x and different valeus n 1 /n 2 are carried out for each bootstrap sample of gauge configurations.The results are quoted from the median and 68% confidence limit of the distribution of all bootstrap samples.Thus, our quoted errors include the correlated statistical and systematic errors arising from x and P z averaging.
Our results for the CS kernel are shown as the black points in Fig. 4. The error bars indicate errors when n 1 /n 2 = 6/7 and 7/8 are excluded from the average.The averages including n 1 /n 2 = 6/7 and 7/8 are depicted as black patches under the data points.
The CS kernel extracted from the GI quasi-TMDWFs calculated in this work is also shown as the blue points and patches, which is consistent with the CG case at smaller b ⊥ .We do not show the CS kernel from the GI quasi-TMDWFs at b ⊥ > 4a because the results are too noisy for comparison as already indicated in Fig. 2. It has been demonstrated in Ref. [55] that after the matching correction, which takes into account of the power corrections at small b ⊥ [54] and the so-called linear renormalon subtraction at large b ⊥ [71], the imaginary part of the CS kernel in the GI case is consistent with zero, so we only take the real part of the final result.
tion [81,82] at the short distances (b ⊥ ≲ 0.4 fm).Beyond this point, the perturbative prediction becomes sensitive to the Landau pole and, thereby, loses reliability.
Although our CS kernel from the GI case loses signal b ⊥ ≳ 0.4 fm, our CG results continue to show very good signals at b ⊥ up to about 1 fm.
For comparisons, we show the most recent lattice QCD calculation (ASWZ24) from GI quasi-TMDWFs in the continuum limit with high statistics [55].Evidently, the results from CG and GI case are consistent with each other, suggesting they fall into the same universality class in the large P z limit under the framework of LaMET [69].Furthermore, our nonperturbative theoretical predictions of the CS kernel are in agreement with the recent phenomenological parameterizations of experimental data, MAP22 [26], ART23 [27], IFY23 [28], and HSO24 (E605) [83,29], which shows a near-linear b ⊥ dependence as proposed in Ref. [84].

Conclusion
We conducted the first lattice QCD calculation of the CS kernel utilizing the recently proposed CG quasi-TMD approach as well as employing unitary domain wall fermion discretization with physical quark masses and a fine lattice spacing.
The CG approach shows significantly lower signal decay, allowing the CS kernel to be determined for extended transverse separations.At the same time, we show that our results are well compatible with the widely used gauge-invariant method.Notably, our results agree with the recent phenomenological parameterizations of experimental data and suggest a near-linear dependence of the CS kernel on large b ⊥ .
This work lays a solid foundation for future research into the CS kernel at larger values of b ⊥ , and advances the QCD computations in the nonperturbative regime of TMD physics.

Figure 1 :
Figure 1: R(t s ; b ⊥ , b z , P z ) as a function of t s for n z = 8 and b ⊥ = 10a.The bands are results from the two-state fits.

Figure 2 :
Figure 2: The real (left) and imaginary (right) parts of the renormalized quasi-TMDWF matrix elements at n z = 8 with b ⊥ = 2a, 6a, 8a for the CG (filled squared symbols) and GI cases (open circled symbols).

Figure 4 :
Figure4: The CS kernel determined from the CG quasi-TMDWFs are shown as the black points and patches, which represent the exclusion and inclusion of momentum pairs n 1 /n 2 = 6/7 and 7/8, respectively.The averages including n 1 /n 2 = 6/7 and 7/8 are depicted as black patches under the data points.The results from GI quasi-TMDWFs calculated in this work are shown as the blue points and patches.For comparison, the CS kernels from recent phenomenological parameterizations of experimental data are shown from MAP22[26], ART23[27], IFY23[28] and HSO24 (E605)[83,29].We also show the perturbative results (N 3 LL) from Ref.[81,82], as well as a recent lattice calculation (ASWZ24) from GI quasi-TMDWFs in the continuum limit with high statistics[55].
program at the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725; the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract DE-AC02-05CH11231 using NERSC award NP-ERCAP0028137.Part of the data analysis is carried out on Swing, a high-performance computing cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory.