Lattice Calculation of the Intrinsic Soft Function and the Collins-Soper Kernel

: We calculate the soft function


Introduction
The transverse-momentum-dependent parton distribution functions (TMDPDFs) [1][2][3], which encode the probability density for 3D parton momenta in hadrons, have been a topic of intense study in modern hadron physics [4,5].The TMDPDFs are universal functions, meaning that they are the same for Drell-Yan (DY) and semi-inclusive deep-inelastic scattering (SIDIS) processes [6], up to at most a sign.Both kinds of experiments have been extensively conducted in the past decades, making up our main knowledge for TMDPDFs [7].The study of TMDPDFs has a long history, including perturbative, phenomenological and non-perturbative determinations, see [8][9][10][11][12] for a selection of recent publications.TMDPDFs can be obtained from experiments by analyzing the final state particles' transverse momenta phenomenologically [6,10,[12][13][14][15][16][17][18][19][20][21][22][23].Such fits always require some non-trivial selection of data, see e.g.Fig. 3 and Tab. 3 in [23] for a recent example.The hard scale Q has to justify the use of perturbation theory and perturbative factorization while, e.g., TMDPDFs are non-perturbative objects depending on q ⊥ and its Fourier conjugate b ⊥ respectively.Although in this case there exists a large amount of data, the resulting error bands in Fig. 8 are large (labeled "ART23").For other TMDs the experimental data situation is much worse, see, e.g., [24].Therefore, combining experimental data with lattice QCD results probably provides the only realistic chance to, e.g., fully determine all eight leading twist TMDs of a nucleon.Such Lattice QCD calculations for TMDs can be grouped in two types.One is based on standard operator product expansion (OPE) plus some additional assumptions to calculate a limited number of Mellin moments of the ratios of TMDPDFs [25][26][27][28].The other follows one of a number of relatively new, more or less equivalent approaches, of which we use the framework of Large Momentum Effective Theory (LaMET) [29,30].In addition to TMDPDF, TMD wave functions (TMDWF) are another important quantity in hadron physics, especially for the description of exclusive reactions.TMDWF provides a description for the partonic structure of hadrons in terms of probability amplitudes.It can be obtained using lattice QCD and LaMET as well [31,32].
LaMET is based on the observation that parton light cone correlations in the rest frame of the hadron, can be obtained from time-independent spatial correlations in the infinite-momentum frame by continuum perturbation theory.At finite but large hadron momenta, LaMET provides a systematic way to determine TMDs via lattice simulations.To do so the universal soft function, which is the focus of this contribution, plays an important role [33].LaMET greatly expands the application of lattice QCD in hadron physics, as reviewed, e.g., in [34].The soft function accounts for non-cancelling soft gluon-radiation at fixed Q ⊥ [1].It consists of a rapidity-independent part called intrinsic soft function S I [35] and a rapidity evolution kernel called Collins-Soper (CS) kernel K [2].
The intrinsic soft function was first introduced in [36] to eliminate the regulator-schemedependence of the quasi TMDPDF/TMDWF.It can be accessed either from heavy quark effective theory [36] or via large-momentum-transfer form factors of light mesons [36].The latter possibility has been explored on the lattice using tree level matching [35,37].The intrinsic soft function was also calculated perturbatively at one-loop order recently in [38].In this work, we will compare the extracted intrinsic soft function using lattice QCD for two different ensembles from the CLS and MILC collaborations.We impose proper normalization and include the one-loop contributions for the first time in a lattice QCD determination.We also apply Fierz rearrangement to suppress higher-twist contaminations [37].
The CS kernel can be obtained from global fits of scattering TMDPDFs data collected primarily for DY and SIDIS processes.It can also be extracted from lattice calculable ratios of either Mellin moments of quasi TMDPDFs/beam functions [39][40][41][42][43] or TMD wave functions (WFs) via a matching procedure.Both tree-level matching [31,35,37] and one loop matching [32] have been explored.In this work, we extract the CS kernel for a CLS ensemble in the framework of LaMET, trying to include one-loop contributions as in [32].Besides, we investigate the pros and cons of extracting the CS kernel from quasi TMDWFs and quasi TMDPDFs.We also compare our results with previous ones, based on experimental [6,12,44,45] and lattice data [31,41,43].
With the CS kernel and intrinsic soft function, a lattice determination of physical TMD-WFs/TMDPDFs based on the factorization Eq.(2.1) and Eq.(2.3) becomes feasible.In [32] the physical TMDWFs are calculated for the first time and in [46] the physical unpolarized TMDPDF of the proton is investigated for the first time on a MILC ensemble.In this work, we discuss TMD-WFs and TMDPDFs as applications of the soft function in TMD physics.We estimate the size of discretization uncertainties by comparing the results for the MILC and CLS ensemble.
The paper is structured as follows: In Sec. 2 we give the theoretical framework of this work.In Sec. 3 we provide the details for the calculation of an intrinsic soft function and in Sec. 4 for the CS kernel.We discuss the application of the soft function for TMDWFs and TMDPDFs in Sec. 5. A summary is given in Sec. 6.

Theoretical framework and lattice setup
The LaMET factorization formula that relates the quasi TMDPDF f to the light cone TMDPDF f reads [47][48][49] where x denotes the longitudinal momentum fraction which is the Fourier conjugate to zP z with z being the offset of the quark-antiquark pair in longitudinal direction and P z the hadron momentum.b ⊥ is the transverse separation that is Fourier conjugate to the transverse momentum Q ⊥ .M is the hadron mass and ζ is a reference rapidity scale that can be chosen at will.The high power corrections are suppressed by the (inverse) hadron momentum and are collected in O(...).Quasi TMDPDFs are also functions of renormalization scale µ and ζ z = (2xP z ) 2 .The hard kernel function H Γ = e h is known at next-to-leading order for TMDPDFs [48] As found in Ref. [38] the matching of TMDWFs is similar to Eq. (2.1) but slightly modified.
Note that following conventions in the literature we use ζ z = (2P z ) 2 for quasi TMDWFs.The hard kernel function H ± = e h ± at next-to-leading order has also a different form [34,38,50] where The subscript ± corresponds to the direction of the Wilson line in quasi TMDWF.
The renormalized quasi TMDWF in momentum space is defined as [51] Ψ± a correlation function which is constructed by inserting a non-local quark-antiquark (q q) current between the vacuum and an external pion state |π⟩.The current is built by connecting the quarkantiquark pair by a staple-shaped Wilson link of length L (or L + |z| for the longer leg) and width b ⊥ where U µ (x; l) ≡ U µ (x, x + ln µ ) and L± ≡ ±max(L, L ∓ z).Fig. 1 depicts how the non-local current is structured: the quark-antiquark pair is shown as black dots connected by the staple-shaped Wilson line shown as thick blue/red lines stretching in +/− direction.If the longitudinal Wilson line points into the positive direction we put use the "+" sign as superscript in Eq. (2.5), otherwise we put "−".In the folowing we take the latter as an example to illustrate our analysis.
< l a t e x i t s h a 1 _ b a s e 6 4 = " A sketch of the non-local quark-antiquark (black dots) current stretching in +/− direction (assuming the momentum to point into the positive z direction).The offset of the two quarks in z-direction is shown as dashed lines of length z.
The determination of TMDWF from the quasi ones requires the intrinsic soft function.The intrinsic soft function was first proposed in LaMET to deal with the divergence related to the emission of soft gluons, which is not cancelled by the real and virtual perturbative corrections.Fortunately according to LaMET it can be isolated and turns into an intrinsic function that can be determined non-perturbatively at small transverse momentum using lattice QCD [34].[36,38] have established a general approach to determine the intrinsic soft function using the quasi TMDWFs Ψ± and a pseudo-scalar light-meson form factor of a transversely-separated product of currents [34,36,38]: P 1 and P 2 denote momenta in opposite directions along the z-axis and are always of equal size in our calculations.The different choices for Γ (=Γ ′ ) project out contributions of different twists, an issue we will address in the next section.The intrinsic soft function then reads [38] S where H (x 1 , x 2 , Γ) is another kernel function known at one loop order [38,50].For Γ = 1 or Γ = γ 5 it reads (2.10) while for Γ = γ ⊥ or Γ = γ ⊥ γ 5 it reads where (2.12) To extract the intrinsic soft function using lattice QCD, a precise determination of the form factor and a well defined quasi TMDWF are necessary.Determing the intrinsic soft function with controlled systematics has developed into a pressing task, in order to expand the range of LaMET applications [35,37,38].
The physical TMDWFs evolve with the rapidity scale ζ satisfying the following renormalization group equation [1,2] which in turn provides the simplest way to determine the CS kernel K(b, µ) from physical TMD data.Physical TMDWFs can be transformed into quasi TMDWFs based on Eq. (2.3).Solving the evolution equation along a constant path of µ allows to fix the CS kernel from the ratio of quasi TMDWFs at different large momenta.The resulting factorization reads [31] which requires a determination of the (renormalized) quasi TMDWFs Ψ± (x, b ⊥ , µ, P z ) on the lattice.Above argument also holds for TMDPDFs, if one simply replaces the (quasi) TMDWFs with the (quasi) TMDPDFs and the accompanying hard kernel function.In fact, the CS kernel is a fundamental nonperturbative function in QCD which describes the interaction of a parton with the QCD vacuum [52].It is believed to be independent of all quantum numbers except for the color representation of the probe.At small b ⊥ , the CS kernel can be reliably determined by perturbative or phenomenological calculations.However at large b ⊥ where the CS kernel becomes non-perturbative, lattice QCD is the only tool that can handle the situation and lattice QCD is essential to relate TMDs at different scales and provides most valuable complementary information compensating lacking experimental data.
With the intrinsic soft function and quasi TMDs at hand, we capture the correct IR physics to all-orders [36,48] and by a perturbative matching the physical TMDs can be obtained.We will illustrate this matching procedure at the end of this paper with two examples, one for a TMDWF and the other for a TMDPDF.
Before diving into the concrete calculations we would like to provide information on the ensembles used throughout this paper in Tab. 1.We use four different ensembles in total.The two CLS ensembles are generated using 2+1 flavor dynamical clover fermions and tree-level Symanzik gauge action.X650 has the same parameters as A654 except for its eightfold larger spatial volume.Note that the light quark and strange quark have the same sea quark mass for these two ensembles.The two MILC ensembles are generated using 2+1+1 flavors of highly improved staggered dynamic quarks [53].They also only differ by their spacial volume which is eightfold larger for a12m130.These ensembles are used in different scenarios: on X650 and a12m310 quasi TMDWFs and form factors are calculated; on A654 and a12m130 quasi TMDPDFs are calculated.In all cases the valence quarks are chosen heavier than the sea quarks for the sake of better signals.The difference due to different spatial size and/or valence/sea masses should be minor [32,37], but will be explicitely investigated in future work.To further improve the signal, hypercubic (HYP) smeared fat links [54] have been used for the staple links in all calculations.In addition, the momentum smearing technique [55] has been used when calculating the quasi TMDPDFs and Coulomb gauge fixed wall source propagators are used when calculating the quasi TMDWFs and form factors.The Table 1.The lattice setups used in this work.X650 and A654 were generated using 2+1 flavors of dynamical clover fermions and tree-level Symanzik gauge action by the CLS collaboration.We remark that for these two ensembles the light quark and strange quark have the same mass in the sea.a12m310 and a12m130 are generated using 2+1+1 flavors of highly improved staggered dynamic quarks [53] (HISQ) by MILC collaboration [56].
last column of the table gives the number of measurements, which is equal to the number of the configurations times the number of different sources used for each configuration.

Intrinsic soft function
As the rapidity independent part of the off-light-cone soft function, the intrinsic soft function S I eliminates the regulator scheme dependence of the quasi TMDPDF/TMDWF.Its determination relies on the calculation of the quasi TMDWF which we present below.

Quasi TMDWF
Bare quasi TMDWF.-FromEq. (2.9) we know that the first piece needed for the intrinsic soft function is the quasi TMDWF.In this section we show how it is determined, taking X650 (a = 0.098 fm) as an example.Similar results have been obtained for a12m130 [31] and a12m310 [32] both with valence pion mass 670 MeV.To obtain the bare quasi-TMDWF in position space on the lattice, one first calculates the following two-point correlation where the interpolators read Ideally one should use Γ = (γ z γ 5 + γ t γ 5 )/2 to eliminate power corrections.However in [31] it was demonstrated that the corrections are at most 5%, such that for simplicity we can just take Γ = γ t γ 5 .In this calculation we have 0 ≤ L ≤ 10a, 0 ≤ z ≤ 10a, 0 ≤ b ⊥ ≤ 7a, 0 ≤ t ≤ 9a and P z = {0, 6, 8, 10, 12} × 2π/(48a) = {0, 1.58, 2.11, 2.64, 3.16} GeV.To ensure that artifacts are small for the considered momenta we examine the dispersion relation in Appendix A for the ensemble X650 and a12m310, on which the soft function will be calculated.We also point out that the previous calculation on a12m130 [31] only considered L = 7a, which should suffice as shall be seen later.We normalize the above non-local two-point function with the corresponding local two-point function and find that the ground-state contribution, which reproduces the continuum definition Eq (2.6), can be obtained by a one-or two-state fit.In [31,32] both fitting Ansätze are explored and a  one-state ansatz was adopted in the end for a better control on the systematic uncertainty in the fits.It is also done here for the same reason.See Fig. 2 for an example for such a fit, where we have defined  Renormalization.-Thebare quasi TMDWFs contain three divergences, the linear divergence originates from the self-energy corrections of the Wilson line, the pinch-pole singularity is caused by the interaction between two legs of the staple-shaped Wilson link and the logarithm divergence is generated by vertices involving Wilson line and light quark.These singularities can be regulated in the way proposed in [51] given by the second line of Eq. (2.5).The square root of the Wilson loop Z E (2L + |z|, b ⊥ , a) renormalizes the former two singularities [57][58][59][60][61] and Z O (1/a, µ) renormalizes the last one [51,62,63].In Fig. 3 we show the Wilson loop calculated on X650 and its extrapolation to large 2L + |z|, where lattice data is either unavailable or too noisy.The extrapolation is feasible because the linear divergence induced by self-energy corrections and gluon exchange effects are exponentially in 2L+|z| [57] and thus can be separated from the rest.The points in the figure denote lattice data and the solid lines are the extrapolated results via one-state fits in the range where we have precise data.We ignore the uncertainties in the extrapolated results as they are negligible compared to other uncertainties, e.g. the statistical uncertainties in the two-point functions.The curves in Fig. 3 actually contain errors, but these are too small to be visible.Arguably the large L limit in Eq. (2.5) can be achieved by looking for a plateau in L. Inspired by the discovery in Ref. [51] the ratio is expected to saturate to a constant at L ≃ 0.7 fm.This does turn out to be the case for our data, see Fig. 4, where a plateau can be identified in the range L ≥ 7a = 0.7 fm for a randomly choosen momentum, b ⊥ and z.In fact, we notice that in this figure the plateau appears already at L ∼ 0.4 fm.As is shown in Appendix C the plateau appears at larger L for smaller z, and L ≃ 0.7 fm is always safe, even at z = 0.  Z O (1/a, µ) can be obtained by taking the ratio of the bare quasi TMDWF calculated at rest on the lattice to the one calculated in the MS scheme where [51] ΨMS The MS scale is set to µ=2 GeV.z 0 and b ⊥0 should be fixed to appropriate values where both discretization artifacts and higher twist contaminations get strongly suppressed, indicated by a Large λ extrapolation.-TheFourier transformation in Eq. (2.5) gets contributions from all real z, including very large ones for which lattice simulations are not possible.For this reason an extrapolation to large z (or equivalently λ) is essential.In Fig. 6 we show the real part of the renormalized quasi TMDWFs at available λ = zP z for a selected momentum P z = 2.64 GeV as an example.We can see that the quasi TMDWFs approach zero at large λ (with larger errors though), which indicates a good convergence when transforming to momentum space.At still larger λ we have to extrapolate.We do so, using the following complex ansatz [61] Ψr− extra (λ) = where m 1 , m 2 , n 1 , n 2 are fit parameters.We perform a joint fit of the tails for b ⊥ in the ranges of λ adjusted for different momenta.In the fits n 1 and n 2 are equal and independent of b ⊥ .m 1 and m 2 are complex valued different for different b ⊥ .λ 0 has been set to a large number, safely larger than the possible correlation length at any finite momentum considered here.A detailed discussion of each term in this ansatz can be found in [46,61].We shift the fit ranges by 1a and reperform the fits.The differences of the central values between the two fits are taken as systematical uncertainties.After λ-extrapolation we Fourier-transform the quasi TMDWFs using The obtained quasi TMDWFs in momentum space are shown in Fig. 7 for a moderate b ⊥ and will be used in the next sections for further calculations.

Pseudo-scalar Meson Form Factor
Extraction of form factor.-Another piece appearing in the definition of the intrinsic soft function is the light pseudo-scalar meson form factor.In this section we calculate this form factor for the two ensembles a12m310 and X650.We choose a12m310 instead of a12m130 based on the practical consideration of computation costs.But we have confirmed that the sea-quark mass effects are marginal, see Appendix D. To allow for large momentum extrapolation we consider the three hadron momenta P z = {4, 5, 6} × 2π/(24a) = {1.72,2.15, 2.58} GeV for a12m310 (a = 0.121 fm) and P z = {6, 8, 10} × 2π/(48a) = {1.58,2.11, 2.64} GeV for X650 (a = 0.098 fm).We have tried including a fourth, higher momentum and found that the difference is negligible due to the larger errors for the fourth momentum.To obtain the bare form factor on the lattice, one needs to calculate a three-point function  Here t seq is the source-sink separation (source at origin) which is set to t seq = {6, 7, 8, 9}a on X650 and t seq = {6, 7, 8, 9, 10}a on a12m310.The interpolators ūΓu and dΓd are inserted at time slice t.It can be shown that after inserting single particle intermediate states and taking the large (imaginary) time limit (0 ≪ t ≪ t seq ), this ratio reproduces the continuum definition Eq. (2.8).In practice this requires a two-state fit of the following ratio .11) in which the ground state contribution gives the bare form factor F (b ⊥ , Γ, P z ).To extract the ground-state contribution one can perform a correlated joint fit for different t seq of the ratio and the two-point function, which share the excitation energy ∆E.This fit is repeated for every P z and b ⊥ .We show an example of such a fit for both X650 and a12m310 at a randomly chosen momentum P z and b ⊥ in Fig. 8.In the fits the first two and last two data points have been discarded due to their strong excited-state contamination.A closer look at the impact of small t seq data set to the fit is given in Appendix E. For the three-point function we have taken the sum µ=1,2 C 3 (Γ = γ µ )+C 3 (Γ = γ 5 γ µ ) to suppress higher-twist effects, which will be discussed in detail in the next paragraph.The extracted bare form factors are renormalized using constants taken from [64] for X650 and for a12m310 we find Z V /Z A = 0.94(1), Z S /Z A = 1.11(3) and Z P /Z A = 0.95 (3).
Fierz rearrangement.-Inspired by Ref. [37] one can Fierz-rearrange the four-quark operators to suppress the higher-twist contaminations.There are a few possibilities to do so.For instance one finds that the combination is dominated by the leading-twist contribution γ µ γ 5 , as the second term in the last line ψc γ µ ψ b ψa γ µ ψ d vanishes for a pion state.In the second line we have used the property that F (γ t ), F (γ z ), F (γ t γ 5 ) and F (γ z γ 5 ) vanish for the pion form factor [38].This verifies the advantage in using such a com- bination of the Dirac structures F (γ ⊥ ) and F (γ ⊥ γ 5 ) on the lattice to identify the leading-twist contribution.
Similarly the combination of F (1) and F (γ 5 ) also gives access to the leading-twist contribution.However, F (1) and F (γ 5 ) have an additional UV divergence [38], leading to a strong momentum dependence, such that this combination is less practical to use, see Fig. 9, where we show the intrinsic soft function obtained from form factors calculated on X650 using tree-level matching (let H(x 1 , x 2 , Γ) = 1 in Eq. (2.9)) at two different momenta from all four channels alone, as well as two combinations of different channels mentioned above.Also shown is the 1-loop perturbative result calculated following [38] using the renormalizaton group equation.The error band is determined in the way described in [32].From the figure we can see that the intrinsic soft functions from different single channels show strong variation, and even carry opposite sign, especially at small b ⊥ , which is however consistent with the observation in [37].When the momentum increases from 2.11 GeV to 2.64 GeV, a better convergence can be seen at larger momentum, confirming the need for large momentum to eliminate power corrections.Another observation is that the intrinsic soft functions from the two combinations show much better consistency, demonstrating that the higher twist effects can be significantly reduced by the Fierz rearrangement.Comparing the intrinsic soft functions from 1−γ 5 and γ ⊥ +γ ⊥ γ 5 , it can be seen that the latter shows only a mild dependence on P z due to the absence of the additional UV divergence [38].Therefore, we use this one in the following calculation.

Results
In Fig. 10 we show the intrinsic soft functions calculated on a12m310 and X650 with treelevel matching and 1-loop matching.Note that the infinite momentum limit is reached only by extrapolation using S I (b ⊥ , µ) = S I (b ⊥ , µ, P z ) + c .The intrinsic soft functions obtained from the combination of F (γ ⊥ ) and F (γ ⊥ γ5) (see Eq. (3.12)).The left panels are from ensemble a12m310 (MILC) while the right panels are from ensemble X650 (CLS).The first row shows results using tree-level matching while the second row is for 1-loop matching.The "−" and "+" sign in the legends indicates the direction of the quasi TMDWF used in the calculation.The data points have been shifted horizontally for better visibility.Only statistical uncertainties are shown.
In all cases the intrinsic soft functions obtained for X650 show stronger P z -dependence than those for a12m310.When going from tree-level matching to 1-loop matching, the intrinsic soft functions increase significantly for both ensembles, approaching the 1-loop perturbative values, especially at small b ⊥ .Based on all these studies, we regard the results from 1-loop matching and γ ⊥ + γ ⊥ γ 5 combination as our final estimates of the intrinsic soft function, summarized in Fig. 11.Generally speaking, the final intrinsic soft function on two ensembles show satisfactory agreement except at small b ⊥ where lattice discretization effects are the most significant.

Collins-Soper Kernel
The CS kernel describes the rapidity evolution of TMDWFs and TMDPDFs.Results containing one-loop contributions were already calculated for a12m130 using quasi TMDWFs in the framework of LaMET in [31] and were revisited in [32] on a12m310.In this section we provide the results for X650 obtained in the same way.Here we use quasi TMDWF in "−" direction and use the 1-loop determination of H [34,38]. usinging the quasi TMDWFs obtained above for different momenta {P z 1 , P z 2 } in Eq. (2.14) we get a momentum-dependent CS kernel.Figure 12.The momentum-dependent and the momentum-independent fits for the CS kernel using the ansatz Eq. (4.1).We select results obtained at b = 0.3fm for ensemble X650 as an example.Only the real parts and statistical uncertainties are shown.
In LaMET in principle both momenta should be large enough to significantly suppress the power corrections.For this reason we choose P z 1 /P z 2 = 3.16 GeV/2.11GeV, 3.16 GeV/2.64GeV.To further extract the leading power contributions, namely to get rid of the finite momentum effects, we fit the momentum-dependent CS kernel to the following theoretically inspired ansatz [31] in the range x ∈ [0.1, 0.9].The intervals beyond this range are discarded as LaMET breaks down there.We show an example of this fit in Fig. 12 for a selected b ⊥ .We point out that the CS kernel calculated in this way is complex and in Fig. 12 only the real part is shown.The imaginary part comes from the matching kernel H, not the quasi TMDWF itself, see [31].The final CS kernel result is shown in Fig. 13 as red points.We take the real part as the central values.The statistical uncertainties are shown as inner error bars and the sum of the statistical and systematical uncertainties are shown as outer error bars.The systematical uncertainties are estimated using the measure From [31] we know that the real part is equivalent to the average of the complex CS kernel calculated for both "±" directions, which at the same time eliminates the imaginary part.
In Fig. 13 we compare the CS kernel obtained in this work with those from other calculations, including the 3-loop perturbative results [44,45], the phenomenological extractions, SV19 [6] and MAP22 [12], and the lattice calculations [31,41,43].The calculation [41] is based on the analysis of the quasi pion beam function with leading order matching kernel.The calculation [31] is same as this work but on the MILC ensemble a12m310.The calculation [43] is based on the analysis of the (first) Mellin moments of the quasi TMDPDF, including one-loop contributions as well.It originally contains four data sets, obtained for pion and proton targets with twist-2 and twist-3 quasi TMDPDF operators.Here we have combined them and shown the results in a single band.In principle the CS kernel can also be obtained from the quasi TMDPDF via Eq.(2.14), by replacing the quasi TMDWF objects with the quasi TMDPDF objects, and replacing the hard kernel function for quasi TMDWF with that for quasi TMDPDF.We have tried this and the results are shown in Fig. 14.In the left panel we show the results obtained for the MILC ensemble a12m130 ( f ) and a12m310 ( Ψ) at a moderate b ⊥ and in the right panel we show the results obtained for CLS ensemble A654 ( f ) and X650 ( Ψ).In the left panel Ψ = ( Ψ− + Ψ+ )/2 so the resulting K(b ⊥ , µ, x, P z 1 , P z 2 ) is pure real while in the right panel only the real parts are shown because the data is lacking for Ψ+ .In the right panel K(b ⊥ , µ, x, P z 1 , P z 2 ) from Ψ− (x, b ⊥ ) with momentum pair 2.64 GeV/2.11GeV shows a poor plateau but this improves when moving to larger momentum pair 3.16 GeV/2.11GeV, consistent with the expectation from LaMET.We want to stress that in this figure we only show the real part for the results obtained from quasi TMDWF while the results obtained from quasi TMDPDF are pure real by construction, which can be seen in Sec.5.2.From this figure we find that for both ensembles the ratios from quasi TMDWF show better plateaus in x while the ratios from quasi TMDPDF decay too fast after x ∼ 0.5, which indicates an early breakdown of LaMET in this range.We remark that the difference we see between the results from quasi TMDWFs and quasi TMDPDFs is unlikely caused by the different pion mass or lattice volume.The common behavior we observed for both MILC and CLS ensembles suggests that this is a generic property.Our tests show that using quasi TMDWF to extract the CS kernel will give better defined results, even though it suffers from systematic uncertainties induced by its imaginary part.To conclude, we will use the results obtained from quasi TMDWF (the squares and circles in Fig. 13) as our final estimate for the CS kernel.

Physical TMDs from the soft function
With the intrinsic soft function and the CS kernel obtained in previous sections, we can extract the physical TMDs on the light cone based on the factorization Eq. (2.1) or Eq.(2.3).We will consider first the TMDWFs in Sec.5.1 and then TMDPDFs in Sec.5.2.Inverting Eq. (2.3) one can obtain the physical TMDWFs.This is done for a12m310 (a = 0.121 fm) and X650 (a = 0.098 fm).The physical TMDWFs are extrapolated to infinite momentum using the ansatz Eq. (3.14) with the intrinsic soft function replaced by the physical TMDWF.The matching was done at rapidity scale ζ = (6 GeV) 2 and renormalization scale µ = 2 GeV.There are various uncertainties contributing to the final physical TMDWFs.They are calculated in the same way as in [32].

TMDWFs in light cone
We show the physical TMDWFs in Fig. 15.The left panel shows the real parts and the right panel shows the imaginary parts.Note that as an example we only show the "−" direction at a small and a large b ⊥ .We remark that we have interpolated the X650 results using a cubic spline to the b ⊥ values analyzed for a12m310 to allow for a direct comparison.The top panels show that there exists visible discrepancies between X650 and a12m310 for the real parts at the two b ⊥ values (in fact, at intermediate b ⊥ values as well, which is not shown here) while for the imaginary part the tension is much reduced.A closer look at these discrepancies requires further investigation of the continuum limit.
Another finding is that the real parts of the amplitudes for both ensembles decrease with b ⊥ , while for the imaginary parts it is the opposite.

TMDPDFs from the light cone
Another application of the soft function is to determine the physical TMDPDFs.This has been done in [46] for a12m130 (a = 0.121 fm), see Tab.1.In this work we show the results obtained for the CLS ensemble A654 (a = 0.098 fm), aiming to understand the lattice discretization effects.The unpolarized bare quasi TMDPDF relevant for this calculations can be found in [46].It should be mentioned that in this work we choose γ t in the bare quasi TMDPDF matrix element, as it approaches γ + at large momentum with smaller operator mixing effects [46].
In the lattice simulations we put 20 sources on each configuration.L is at most 10a and b ⊥ /a is at most 7a.Both can be positive or negative.We average the two directions because they are equivalent.z is at most 15a and can also point into both directions.The real part of the threepoint function is symmetric with respect to z = 0 so it is averaged again.However the imaginary part is anti-symmetric so we need the multiply the negative part by −1 before averaging.In this way, we have 160 measurements per configuration.To control the excited states contamination, we calculate correlators at four different source-sink separations t seq /a = {5, 6, 7, 8} and extract the ground state contribution through a joint fit for all separations.The fits are always performed in the range t ∈ [1, t seq − 1].To examine the dependence on momentum, we consider three nucleon momenta P z ∈ {3, 4, 5} × 2π/(24a) = {1.58,2.11, 2.64} GeV.
The bare quasi TMDPDF is renormalized in the same way as the quasi TMDWF using Wilson loop and Z O [51].The difference is that Z O is now determined from the TMDPDF calculated in MS scheme hMS Γ (z 0 , b ⊥0 , µ), which turns out to be the same as the TMDWF in MS scheme at zero momentum.On A654 the best plateau in b 0 appears at z 0 = 0, starting from b ⊥0 = 1a = 0.1 fm and we choose z 0 = 0, b ⊥0 = 1a.With proper renormalization group evolution (RGE) as done in [46], we find Z O = 1.406 (1).As for the large L limit, following the spirit of Ref. [51] and as verified above, we take the data at L = 7a = 0.7 fm at which the renormalized quasi TMDPDF has saturated to a constant.
For the λ extrapolation, which uses the ansatz Eq. (3.8), we first fit the tails in the range λ ∈ [9a, 12a]P z at all momenta and estimate the central values and statistical uncertainties from this fit.Then we repeat the fit in the range λ ∈ [7a, 10a]P z .The differences of the central values between the two fits are taken as the systematical uncertainties of the extrapolation.After λextrapolation we Fourier-transform the matrix elements obtained in position space to momentum space.We remark that the imaginary part vanishes due to its antisymmetry with respect to z = 0.
We now calculate the physical TMDPDF by applying the matching Eq. (2.1) using the scale ζ = 4 GeV 2 , followed by an extrapolation to infinite momentum.As in [46], the RGE for the hard kernel Eq (2.2) is solved.In Appendix F we illustrate the evolution of TMDWF and TMDPDF with respect to the renormalization scale µ and rapidity scale ζ.
In Fig. 16 we show the physical TMDPDFs obtained on A654 and compare them to those obtained on a12m130 [46].The A654 results are again interpolated to the b ⊥ values attainable on a12m130.The endpoint regions x < 0.2 and x > 0.8 where LaMET breaks down have been grey shaded.Due to renormalization group resummation, the shaded regions are broader than in the case of the TMDWF.Overall speaking, the TMDPDFs on different ensembles have similar shape at the two b ⊥ values considered here, in fact at other b ⊥ as well, which we do not show; but given the uncertainties at b ⊥ = 0.12 fm and the observed discrepancies, further study of in particular discretization effects is needed to obtain reliable error estimates.

Conclusion
We calculate the intrinsic soft function and the CS kernel on two different lattice ensembles.In this updated result of the soft function we have included the one-loop contributions and used proper normalization of light meson form factors.We have also suppressed higher-twist contaminations by Fierz rearrangements.We find that the intrinsic soft function obtained on two ensembles are similar except at small b ⊥ where lattice discretization effects are probably significant.
We have also tried to extract the CS kernel from TMDWFs including the new X650 ensemble in addition to those studied in [31].We also extract the CS kernel from quasi TMDPDFs.It turns out that using the former method is a better choice.We provide a comparison of the CS kernel obtained in this work and in other studies.We find that the CS kernel calculated on X650 shows good agreement with literature, particularly [31].Using the soft function we determine the physical TMDWFs for pion and physical TMDPDFs for proton from the corresponding quasi TMD objects renormalized using the method proposed in [51], also on two ensembles.The good agreement observed on different ensembles for the TMDWFs and TMDPDFs verifies the applicability of calculating light cone quantities from lattice simulated quasi objects using TMD factorization via the soft function.From our findings we conclude that a determination of the soft function with better controlled precision and a systematic investigation of discretization effects is needed and possible.We leave this for future work.sufficient.Comparing the top panels and the bottom panels, one finds that for larger z, it is still safe to include more data points from smaller L in the fit.

D Impact of sea quark mass to the determination of S I
To examine the impact of sea quark mass on the determination of S I , we calculate the form factor and quasi TMDWF at a selected momentum P z =1.72 GeV on a12m130 and a12m310.For brevity, in the calculation of the form factor we consider only one single t seq = 8a and take the value at t seq /2 as an estimate for the t seq → ∞ result, instead of performing a real fit.The results are shown in Fig. 20.It can be seen that for both quantities results from both ensembles give consistent results.The tiny difference is much smaller than the uncertainties from other aspects, e.g., large-λ extrapolation and the matching procedure, and thus can be safely ignored.E Impact of small t seq to the fit of three-point correlation function To see how much the small t seq data set affects the extrapolation to t seq → ∞, we perform a fit of the same data set used in the right panel of Fig. 8 but exclude t seq = 6a.The results are shown in Fig. 21.Comparing to the right panel of Fig. 8 it can be seen that only the error in F (b ⊥ , µ) grows slightly while the change in the mean value lies within the statistical error.This confirms that the excited-state contamination is under control after the extrapolation.If we include the change coming from different data sets used in the extrapolation as systematic uncertainty, we can see it is of O(1%), similar to the statistical uncertainty.Such change will be much smaller for the MILC ensemble where more t seq are available.We stress that its influence is very limited in the sense of physics, when compared to, e.g., the change of matching from tree level to one-loop level, see Fig.

F Evolution of TMDWF and TMDPDF with ζ and µ
In this section we illustrate the evolution of TMDWF and TMDPDF with respect to the rapidity scale ζ and the renormalization scale µ.For TMDWF we take a MILC ensemble as an example.We consider three renormalization scales µ = 1.5, 2.0, 2.5 GeV and three rapidity scales ζ = 25, 36, 49 GeV 2 .The results are shown in Fig. 22.For TMDPDF we take a CLS ensemble as an example.The renormalization scales are the same as above and the rapidity scales are ζ = 3, 4, 5 GeV 2 .The results are shown in Fig. 23.It can be seen that in all cases TMDWF and TMDPDF have very mild dependence on both scales in the chosen range considered here.
) containing mainly three parts: the bare quasi TMDWF Ψ± (L, z, b ⊥ , P z ) in position space, the Wilson loop Z E of length 2L + |z| and width b ⊥ , and a quark Wilson line vertex renormalization factor Z O .The latter two are part of the renormalization procedure which we will return to in the next section.a denotes the lattice spacing.The bare quasi TMDWF in position space reads

Figure 2 .
Figure 2. Example for a one-state fit in t with momentum P z = 2.64 GeV, b ⊥ = 0.3 fm, L = 0.7 fm and z = 0.3 fm.
) for simplicity.The figure shows that the one-state Ansatz does capture the correct behavior of the lattice data.More examples are given in Appendix B.

Figure 4 .
Figure 4. Fit to a constant at large L for momentum P z = 2.64 GeV, b ⊥ = 0.2 fm and z = 0.4 fm.

Figure 7 .
Figure 7.The real (left) and imaginary (right) part of the quasi TMDWFs obtained for X650 in momentum space with different P z for a selected b ⊥ .

Figure 8 .
Figure 8. Extraction of the bare form factor via a joint fit for the ratio Eq. (3.11) and the local two-point function Eq. (3.1).The data points are the lattice data of the ratio and the colorful bands denote the fit results.The ground-state contributions (namely the bare form factors) are shown as grey bands.The left panel is for ensemble a12m310 (MILC) and the right panel is for ensemble X650 (CLS).

Figure 9 .
Figure 9.The intrinsic soft function calculated on X650 using tree-level matching from different channels and two combinations of them.The left panel shows results obtained at 2.11 GeV while the right panel shows results obtained at 2.64 GeV.The data points have been shifted horizontally for better visibility.Only statistical uncertainties are shown.
Figure 10.The intrinsic soft functions obtained from the combination of F (γ ⊥ ) and F (γ ⊥ γ5) (see Eq. (3.12)).The left panels are from ensemble a12m310 (MILC) while the right panels are from ensemble X650 (CLS).The first row shows results using tree-level matching while the second row is for 1-loop matching.The "−" and "+" sign in the legends indicates the direction of the quasi TMDWF used in the calculation.The data points have been shifted horizontally for better visibility.Only statistical uncertainties are shown.

Figure 11 .
Figure 11.Final results for the intrinsic soft functions obtained from the CLS and MILC ensembles.S lat,1 loop± corresponds to the lattice results extracted by Ψ± .

Figure 14 .
Figure 14.The momentum-dependent CS kernel obtained from the quasi TMDWF and quasi TMDPDF on MILC ensembles (left) and CLS ensemble (right).In both panels the results are obtained from different ensembles: a12m130 for Ψ and f in the left panel, as well as X650 for Ψ and A654 for f in the right panel.

Figure 15 .
Figure15.The panel shows the real parts of the physical TMDWFs from a12m310 and X650 while the right panel show the imaginary parts.They are determined at ζ = (6 GeV) 2 , µ = 2 GeV and in the infinite P z limit.Only Ψ − is considered as an example.

xfFigure 16 .
Figure 16.A comparison of the physical TMDPDFs obtained on CLS ensemble A654 and MILC ensemble a12m130 at b ⊥ =0.12 fm (left) and 0.48 fm (right).On A654 we have interpolated TMDPDFs using a cubic spline to the b ⊥ values used for a12m130 for a direct comparison.

z 8 zzzFigure 19 .
Figure 19.More examples with {z, b } for large L fit at P z = 2.64 GeV.

Figure 20 .bFigure 21 .
Figure 20.Form factor (left) and quasi TMDWF (right) calculated at Pz=1.72 GeV on a12m130 and a12m310.In the calculation of the form factor we take the value at t = tseq/2 = 4a as an estimate for the ground state contribution.In the right panel we choose a moderate b ⊥ =0.36 fm.