Determination of the Collins-Soper kernel from Lattice QCD

This work presents a determination of the quark Collins-Soper kernel, which relates transverse-momentum-dependent parton distributions (TMDs) at different rapidity scales, using lattice quantum chromodynamics (QCD). This is the first such determination with systematic control of quark mass, operator mixing, and discretization effects. Next-to-next-to-leading logarithmic matching is used to match lattice-calculable distributions to the corresponding TMDs. The continuum-extrapolated lattice QCD results are consistent with several recent phenomenological parameterizations of the Collins-Soper kernel and are precise enough to disfavor other parameterizations.

In this context, a quantity of particular importance is the Collins-Soper (CS) kernel: a fundamental nonperturbative function that appears as the universal rapidity evolution kernel for TMDs, which can be considered to characterise the QCD vacuum [22][23][24].The CS kernel is not only a fundamental proton and nuclear structure observable of importance in its own right, it is also needed to compare TMDs measured at different scales and is required as input into measurements of electroweak observables including the W-boson mass [25] and in various nuclear structure studies [20].
Phenomenological extractions of the CS kernel from global fits of experimental data from Drell-Yan and Semi-Inclusive Deep-Inelastic Scattering (SIDIS) processes, however, are largely unconstrained in the nonperturbative region, specifically for kinematics with small transverse momentum scale q T < ∼ 0.3 GeV.First constraints of the CS kernel from lattice QCD calculations [26][27][28][29][30][31][32][33][34][35] demonstrate the potential of this approach to provide first-principles information with sufficient precision to distinguish between different phenomenological models in this regime.Nevertheless, key systematic uncertainties, in particular discretization effects, remain to be controlled in all such calculations to date.
This work presents a new lattice QCD determination of the CS kernel which includes systematic control of quark mass, operator renormalization, and discretization effects, and uses next-to-next-to-leading logarithmic matching to TMDs from the corresponding lattice-calculable distributions.This allows a parameterization of the CS kernel to be constrained entirely by first-principles calculations for the first time.
The Collins-Soper Kernel: The transverse momen-tum of a parton of flavor i in a given hadron state is encoded in the TMDs f TMD i (x, b T , µ, ζ), which are functions of the longitudinal momentum fraction x carried by the parton, the transverse displacement b T (the Fourier conjugate of q T ), the virtuality scale µ, and the hadron momentum through the rapidity scale ζ.Unlike the µevolution of the TMDs, which is perturbative for perturbative scales µ and ζ, the ζ evolution of TMDs is nonperturbative and is encoded in the CS kernel [22,23]: The quark CS kernel γ q (b T , µ) is independent of flavor.
Lattice QCD calculation: Constraints on the quark CS kernel are extracted from lattice QCD calculations on 3 ensembles of gauge fields produced by the MILC collaboration [42] with 2+1+1 dynamical quark flavors, the one-loop Symanzik improved gauge action [43][44][45][46], and the highly improved staggered quark action with sea quark masses tuned to reproduce the physical pion mass [47][48][49].The calculations are performed as detailed in Ref. [35], which presented results on one ensemble of lattice gauge fields with four-volume L 3 × T = (48a) 3 × 64a with a = 0.12 fm.This work adds calculations on an additional two ensembles of gauge fields with four-volumes L 3 × T = (32a) 3 × 48a with a = 0.15 fm and L 3 × T = (64a) 3 × 96a with a = 0.09 fm, enabling systematic investigation of discretization effects for the first time.A summary of the computation and analysis are included below, with further details and figures provided in the Supplementary Material.
The numerical calculation proceeds as detailed in Ref. [35]: Computation of bare quasi-TMD WF ratios: Bare quasi-TMD WFs φΓ (b T , b z , P z , ℓ, a) are extracted from bootstrap-level fits to the Euclidean time-dependence of two-point correlation functions both with and without staple-shaped operators.Pion states are created with momentum-smeared interpolating fields [70], and the tree-level O(a)-improved Wilson clover fermion action [71][72][73] is used for propagator computation, with clover term coefficient c sw = 1.0.Hopping parameters are κ ∈ {0.12575, 0.12547, 0.1252} for calculations on the ensembles with a ∈ {0.15, 0.12, 0.09} fm respectively, and field configurations are treated with Wilson flow to flow-  L n z .For operators with staple extent ℓ/a, all geometries with −ℓ/a ≤ b z ≤ ℓ/a and 0 ≤ bT /a ≤ 7 along nT ∈ {±x, ±ŷ} are computed, for all of the 16 Dirac structures Γ.The number of configurations used for each measurement is denoted N cfg ; on each configuration, sources on a 2 4 grid bisecting the lattice along each dimension are used.time t = 1.0 [74] and gauge-fixed to Landau gauge1 before measurements are made.
Calculations are performed on each ensemble for the operator choices (defined by Γ and the staple geometry specified by ℓ, b z , b T ), choices of momenta P z , and the numbers of configurations specified in Table I.For the geometries with odd ℓ/a where no b z = 0 matrix elements are available (as used in the denominator of Eq. (3) to form ratios), the average of those with b z /a = ±1 are used.
The real parts of the MS-renormalized quasi-TMD WF ratios, W MS Γ (bT , µ, b z , P z , ℓ, a), computed on each ensemble.Panels show results on ensembles with a ∈ {0.15, 0.12, 0.09} fm from top to bottom, as functions of b z , for Γ = γ4γ5, P z = 1.3 GeV, and similar though not identical bT values as indicated in each panel.
Fourier transformation: After renormalization and multiplication by N Γ (P z ), a Discrete Fourier Transform (DFT) is used to realize the Fourier transforms in the numerator and denominator of Eq. ( 2), where the quasi-TMD WF ratios for each P z are first averaged over ±b z and the relevant values of ℓ(P z ).This yields x-space renormalized quasi-TMD WF ratios W MS Γ (b T , µ, x, P z , a).Because results are computed for a finite range of b z , the DFT is effectively truncated to a finite range.The effects of this truncation are studied by comparing results using subsets of the data with b z < b z max and varying b z max , as well as by comparing results where an analytical model used to extrapolate to b z > b z max .Truncation effects are found in this way to lead to negligible systematic uncertainties, as detailed for the a = 0.12 fm ensemble in Ref. [35].
Perturbative matching: The perturbative matching correction δγ MS q (µ, x, P z 1 , P z 2 ) appearing in Eq. ( 2) is taken to be the 'b T -unexpanded resummed next-to-next-to-leading order' (uNNLL) correction detailed in Ref. [35].It was found in the analysis of Ref. [35], for the same a = 0.12 fm ensemble also studied here, that this choice reduces the effect of b T -dependent power corrections and offers the best convergence compared with other currently-available matching prescriptions, i.e., fixed-order [63,[78][79][80] and resummed [59,64,81] corrections up to next-to-next-toleading order.
Extraction of the CS kernel: Each choice of x, P z , and a defines an estimator for the CS kernel (corresponding to the right-hand side of Eq. ( 2) with the Fourier Transform implemented via DFT, neglecting the p.c. term, and without the limit a → 0 being taken): In principle, it would be desirable to exploit the multiple lattice ensembles available in this calculation to perform a continuum extrapolation of the CS kernel estimator for individual b T values; this would require matched bT across the ensembles.Alternatively, one might aim to disentangle power corrections and discretization effects by fitting all results for Re γMS Γ (b T , µ, x, P z 1 , P z 2 , a) to a parameterization of the CS kernel plus P z , a, and b T -dependent terms 2 .In practice, estimators for different {P z 1 , P z 2 } are largely consistent, and as such the data is insufficient to constrain momentum-dependent power corrections.
Instead, the CS kernel on each ensemble is determined as a bootstrap-level weighted average of Re γMS Γ (b T , µ, x, P z 1 , P z 2 , a) over Γ ∈ {γ 4 γ 5 , γ 3 γ 5 }, all available combinations of {P z 1 , P z 2 }, and x ∈ [0.3, 0.7], with weights proportional to the inverse variance, just as done in Ref. [35].These P z -, Γ-, and x-averaged CS kernel constraints, denoted γ MS q (b T , µ, a), should agree with the CS kernel up to discretization effects.The results (including a fit to a parameterization of the CS kernel and the leading a/b T discretization effects, as discussed in the next section), are shown in Fig. 2.
Additionally, an analogous analysis of Im γMS Γ (b T , µ, x, P z 1 , P z 2 , a) can be performed.As the CS kernel is purely real, significant deviation of the resulting numerical results from zero would provide an indication of systematic uncertainties beyond those that are quantified in this calculation.This analysis is presented in the Supplementary Material.Including the uNNLL matching of Ref. [35] and recent devlopments [82] accounting for a linear infrared renormalon in the imaginary part of the matching coefficient for the quasi-TMD WF, there is no evidence in the numerical data for significant additional unconstrained systematic uncertainties.

Parameterization:
The Lattice QCD constraints on the CS kernel are fit to the parameterization of Ref. [40], with the addition of terms accounting for lattice discretization effects proportional to a/b T , a 2 /b 2 T : where 3 D res is the resummed leading power expression for the CS kernel computed in the operator product expansion, evolved to scale µ, and the parameterization of the remaining nonperturbative piece is 2 Specifically such corrections would be proportional to terms such as a/b T , a 2 (P z 1 ) 2 − a 2 (P z 2 ) , . . . 3 Dres is given explicitly in the Supplementary Material.
[Upper] Averaged CS kernel estimators computed on each ensemble, including a fit to a parameterization of the CS kernel plus O(a/bT ) discretization effects, as described in the text.The colored dashed curves correspond to γ param q (bT , µ, a), with the best-fit values of (B NP , c0, c1, k1, k2) as described in the text, at each corresponding value of a, while the solid black curve shows the result at a = 0. [Lower] Lattice QCD constraints on the CS kernel, with O(a/bT ) artefacts subtracted as defined in the text, and the best-fit parameterization of the CS kernel fit to the lattice results shown as a solid black curve with the 1σ uncertainty band shown as a shaded red region.
The expression of Eq. ( 6) is thus a three-parameter (B NP , c 0 , c 1 ) parameterization of the CS kernel, with an additional two parameters (k 1 , k 2 ) modeling lattice discretization effects.
quantify the relative goodness-of-fit for models including different parameter subsets.The minimum AIC model is found to be (c 0 , k 1 ) with c 1 = k 2 = 0 and B NP = 2 GeV.The corresponding fit results are c 0 = 0.032 (12), k 1 = 0.22 (8), (9) with a χ 2 /dof = 0.39.These fit results, and the resulting parameterization of the CS kernel, are shown in Fig. 2. Overall fit quality is illustrated through the comparison of γ param q (b T , µ, a = 0) with best-fit values for (B NP , c 0 , c 1 , k 1 , k 2 ) with the lattice QCD results where discretization effects have been subtracted, i.e., γ MS q (b T , µ) ≡ γ MS q (b T , µ, a) − k 1 (a/b T ) using the best-fit results for k 1 .
These continuum-limit results are compared with phenomenological parameterizations of experimental data in Fig. 3.In particular, the parameterization used in Ref. [37] corresponds to the AIC-preferred parameterization used here and leads to a consistent result c SV19 0 = 0.043 (11) with B SV19 NP = 1.9(2)GeV.The global fits performed in Ref. [40] also give a consistent result, c ART23 0 = 0.037 (6), though in that work c 1 is also included as a fit parameter.
Fits to other parameter subsets (c 0 , k 2 ) and (c 0 , k 1 , k 2 ) give consistent results for c 0 at 1σ with uncertainties that differ by < ∼ 10%.The magnitudes of k 1 and k 2 range from 0.1 -0.3 in all cases, which suggests that the size of discretization effects is consistent with naive dimensional analysis.Fits including B NP or c 1 as free parameters give consistent results for c 0 with larger uncertainties.
Other parameterizations for the nonperturbative function D NP (b) have been used in fits to experimental data [36,87], for example the BLNY parameterization D BLNY NP (b) = g 2 b 2 with free parameters g 2 and B NP (which enters D res ).Fits to this parameterization with B NP = 1.5 GeV lead to the result g 2 = 0.085 (26) with compara-ble goodness-of-fit, χ 2 /dof = 0.58, to the fits using the parametrization of Eq. (7) described above.This is consistent with the phenomenological fit results of Ref. [41], which use the same value of B NP and find g 2 = 0.053 (24).These lattice QCD constraints on the CS kernel are therefore not sufficient to establish a clear preference between functional forms for D NP ; however they do provide a significant preference for the recent fit results from Refs.[37,[39][40][41] in comparison with Ref. [38] and especially with older BLNY fit results [36] at large b T .
Summary: This work presents the first lattice QCD calculation of the CS kernel with systematic control of quark mass, operator renormalization, and discretization effects.The results are used to constrain a 'pure-theory' parameterization of the CS kernel through a direct fit to lattice QCD results for the first time.These lattice QCD results for the CS kernel are consistent with the most recent phenomenological results.This opens the door for future first-principles QCD predictions of the CS kernel beyond the region constrained by current experiments, as well as joint fits to experimental data and lattice QCD results.As more precise lattice QCD results are achieved at larger values of b T in future calculations, this promises to be increasingly valuable.PES

SUPPLEMENTARY MATERIAL
This Supplementary Material collects additional results of the intermediate stages of analysis in the numerical calculation of the CS kernel.

A. Two-point correlation functions and effective energies
For each of the ensembles detailed in Table I, two-point correlation functions both with and without staple-shaped Wilson line operators are constructed, and used to extract the bare quasi-TMD WF ratios ratios 3), as described in the text.Pion interpolating operators are defined with momentum-smeared quark fields constructed with a Gaussian momentumsmearing kernel F K with K = ±P/2, realized with n smear = 50 smearing steps and width ε = 0.2 [70]: The two-point correlation functions are constructed as where labels a denoting the discretization scales of different ensembles are suppressed, and P = P z ẑ, t = y 4 .
Energies E π (P z ẑ) on each ensemble, for each considered value of momentum P z , are extracted from effective energies constructed using both C π 2pt and the most statistically precise C Γ 2pt for each choice of P z , as detailed in Ref. [35]: Effective energies on each ensemble are shown in Fig. 4, and the extracted values of E π (P z ẑ), compared with the continuum dispersion relation 2 , are shown in Fig. 5. Deviations between E π (P z ẑ) and the continuum dispersion relation are observed to increase approximately linearly with aP z , as expected since tree-level O(a) improvement only partially mitigates O(a) discretization effects.The deviations from the continuum dispersion relation on each ensemble are < ∼ 19% with a = 0.15 fm, < ∼ 15% with a = 0.12 fm, and < ∼ 10% with a = 0.09 fm.
Effective energies (Eq.( 12)), constructed using C π 2pt (squares, offset slightly on the horizontal axis) and the most statistically precise C Γ 2pt (triangles) for each choice of P z on each ensemble.The gray bands show fit results obtained as detailed in Ref. [35] as weighted averages over fits with different fit ranges and numbers of excited states, while the colored bands show the single highest-weight fit included in each average.

B. Bare and renormalized quasi-TMD WF ratios
Examples of bare quasi-TMD WF ratios W The corresponding results for the a = 0.12 fm ensemble are presented in Ref. [35].

Results
for the CS kernel estimator γMS Γ (b T , µ, x, P z 1 , P z 2 , a) (Eq.( 5)) averaged over x ∈ [0.3, 0.7] and Γ ∈ {γ 4 γ 5 , γ 3 γ 5 } but shown separately for each momentum pair P z 1 , P z 2 are presented for all three ensembles in Figs.37-39.Analogous results for the CS kernel estimator averaged over x and P z 1 , P z 2 but shown separately for each Γ are also presented in Figs.[37][38][39] The results shown for the a = 0.12 fm ensemble are reproduced from Ref. [35] to enable comparison with the results from the other two ensembles.

D. Imaginary parts of CS kernel estimators
Previous analyses of lattice QCD results for the CS kernel in Refs.[32,35] yielded nonzero values of the imaginary part.Without explanation, this would indicate the presence of systematic uncertainties unaccounted for in the calculations.However, these observations can be understood in the context of a recent analysis [82] which reveals that the imaginary part of the matching coefficient for the quasi-TMD WF includes a linear infrared renormalon.To mitigate the nonconvergence of the perturbative series, it was proposed in Ref. [82] to subtract the leading asymptotic series from the matching coefficient, i.e., a leading-renormalon subtraction (LRR) scheme, with the uncertainties in the perturbative series absorbed into a linear power correction.
The matching correction to the CS kernel in Eq. ( 5) can be expressed as where x = 1 − x, and the matching coefficient H is the product of two matching kernels for the 'heavy-light' current, which are defined in Eqs.(C2-C8) in Ref. [35].Here the b T -dependence in δγ MS q is suppressed as it exists only in the b T -unexpanded matching coefficient in Eq. (C43) of Ref. [35].
The asymptotic series in C is given by from Eqs. (2.48-52) of Ref. [82], where β 0 is the lowest order β-function, and N m = 0.552 for n f = 4. Therefore, the LRR subtracted matching coefficient is simply which is shown to converge quickly in the perturbative order [82].For the analysis in this work, C is thus replaced by C LRR , including the b T -unexpanded uNNLL matching coefficient, defined in Eq. ( 27) and Eq.(C43) of Ref. [35], to obtain the final quoted results [35].
Since the linear renormalon is a feature in the TMD factorization of the quasi-TMD WF [82], the LRR is only expected to improve the perturbative convergence in the large b T , or large P z b T , region.When P z b T < ∼ 1, the b T -dependent uNNLL matching kernel takes care of the perturbative power corrections, and has been demonstrated to yield an imaginary part of the CS kernel which is consistent with zero [35].As shown in Fig. 40, in the numerical calculation of this work results with both LRR and uNNLL matching indeed display the expected behavior in the relevant parameter regions.Thus, the results do not imply additional uncontrolled systematic uncertainties beyond those discussed in the main text.

E. Perturbative part of the CS kernel
The perturbative part of the CS kernel D res in Eq. ( 6) is given by where µ b * = 2e −γ E /b * , Γ cusp is the cusp anomalous dimension known to four loops [100][101][102][103][104][105] and approximated at five-loop order [106], and d is the non-cusp part of the rapidity anomalous dimension known to four loops [85,86,104,105].
Both anomalous dimensions are perturbative series in the strong coupling α s , where a s = α s /(4π).For N 3 LL resummation, the required anomalous dimensions are and Additionally, the QCD β-function is defined as with The kernel K in Eq. ( 17) at where r = α s (µ)/α s (µ 0 ).The four-loop running coupling is needed for N 3 LL resummation: In the numerical analysis, the initial condition of α s is set to be α s (µ 0 = 2 GeV) = 0.293 and n f = 4.

5 FIG. 5 .
FIG. 5. Extracted energies Eπ(P z ẑ) on each ensemble as functions of P z .The energy determined from the continuum dispersion relation is shown as a grey dashed line.

Γ
(b T , b z , P z , ℓ, a) (Eq.(3)) for all Γ and all relevant ℓ and b z are shown for representative values of b T and P z for all ensembles in Figs.6-8.Renormalized quasi-TMD WF ratios W MS Γ (b T , µ, b z , P z , ℓ, a) for the Γ ∈ {γ 4 γ 5 , γ 3 γ 5 } used to constrain the CS kernel are shown for all b T , b z , P z , and ℓ for the a = 0.15 fm and a = 0.09 fm ensembles in Figs.9-36.
is supported in part by Simons Foundation grant 994314 (Simons Collaboration on Confinement and QCD Strings).YZ is also supported in part by the 2023 Physical Sciences and Engineering (PSE) Early Investigator Named Award program at Argonne National Laboratory.This manuscript has been authored by Fermi Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S. Department of Energy, Office of Science, Office of High Energy Physics.This research used resources of the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231, the Extreme Science and Engineering Discovery Environment (XSEDE) Bridges-2 at the Pittsburgh Supercomputing Center (PSC) through allocation TG-PHY200036, which is supported by National Science Foundation grant number ACI-1548562, facilities of the USQCD Collaboration, which are funded by the Office of Science of the U.S. Department of Energy.We gratefully acknowledge the computing resources provided on Bebop, a high-performance computing cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory, and the computing resources at the MIT SuperCloud and arXiv:1812.11818[hep-ph].[107] A. H. Hoang, D. W. Kolodrubetz, V. Mateu, and I. W. Stewart, Phys.Rev. D 91, 094017 (2015), arXiv:1411.6633[hep-ph].