Nonperturbative Determination of Collins-Soper Kernel from Quasi Transverse-Momentum Dependent Wave Functions

In the framework of large-momentum effective theory at one-loop matching accuracy, we perform a lattice calculation of the Collins-Soper kernel which governs the rapidity evolution of transverse-momentum-dependent (TMD) distributions. We first obtain the quasi TMD wave functions at three different meson momenta on a lattice with valence clover quarks on a dynamical HISQ sea and lattice spacing $a=0.12$~fm from MILC, and renormalize the pertinent linear divergences using Wilson loops. Through one-loop matching to the light-cone wave functions, we determine the Collins-Soper kernel with transverse separation up to 0.6~fm. We study the systematic uncertainties from operator mixing and scale dependence, as well as the impact from higher power corrections. Our results potentially allow for a determination of the soft function and other transverse-momentum dependent quantities at one-loop accuracy.

In the framework of large-momentum effective theory at one-loop matching accuracy, we perform a lattice calculation of the Collins-Soper kernel which governs the rapidity evolution of transversemomentum-dependent (TMD) distributions. We first obtain the quasi TMD wave functions at three different meson momenta on a lattice with valence clover quarks on a dynamical HISQ sea and lattice spacing a = 0.12 fm from MILC, and renormalize the pertinent linear divergences using Wilson loops. Through one-loop matching to the light-cone wave functions, we determine the Collins-Soper kernel with transverse separation up to 0.6 fm. We study the systematic uncertainties from operator mixing and scale dependence, as well as the impact from higher power corrections. Our results potentially allow for a determination of the soft function and other transverse-momentum dependent quantities at one-loop accuracy.

I. INTRODUCTION
Understanding the internal, three-dimensional structure of hadrons, such as the proton, is an important goal in nuclear and particle physics. In this regard, the transverse momentum-dependent (TMD) parton distribution functions (TMDPDFs) [1,2] play an important role as they characterize their intrinsic transverse partonic structure. These distributions are also essential ingredients in the description of multi-scale and noninclusive processes, such as Drell-Yan production of electroweak gauge bosons or Higgs bosons or semi-inclusive In contrast to the TMDPDFs that encode the probability density of parton momenta in hadrons, the transverse momentum-dependent wave functions (TMDWFs) offer a probability amplitude description of the partonic structure of hadrons, from which one can potentially calculate various quark/gluon distributions. In the QCD factorization involving transverse momentum, they are the most important ingredients to predict physical observables in exclusive processes, for instance, weak decays of heavy B meson [8,9] which are valuable to extract the CKM matrix element, and to probe new physics beyond the standard model. However, due to the lack of knowledge of TMDWFs, the one-dimensional lightcone distribution amplitudes (LCDAs) are used instead in most analyses of B decays [8][9][10], resulting in uncontrollable errors. The unprecedented precision of experimental measurements of B decays [11] urgently requires a reliable theoretical knowledge of TMDWFs.
A common feature of TMDPDFs and TMDWFs is that they depend both on the longitudinal momentum fraction x and on the transverse spatial separation of partons. Considerable theoretical efforts have been devoted in recent years to determine these quantities by fitting the pertinent experimental data [12][13][14][15][16][17][18][19][20], which, however, is limited by the imprecise knowledge of the nonperturbative behaviour of TMDPDFs and TMDWFs. Thus, it is highly desirable to develop a method to calculate them from first-principle approaches such as lattice QCD.
This has been realized in the framework of large momentum effective theory [21,22], which offers a systematic way to calculate light-cone correlations by simulating time-independent Euclidean correlations on the lattice. Significant progress has been made in calculating various parton quantities from LaMET. For recent reviews, see Ref. [23,24].
A very important result of LaMET development is that the TMDPDFs and TMDWFs can be calculated through the Euclidean quasi-TMDPDFs and quasi-TMDWFs, as well as a universal soft function (factor) [25][26][27][28]. In Ref. [26], it has been suggested that the form factor of a bi-local four-quark operator calculable on the lattice, can be factorized into quasi-TMDWFs, a universal soft (function) factor and the matching kernel through QCD factorization at large momentum transfer, allowing for the first time calculating the universal soft function on lattice. Thus the light-cone TMD parton distributions and wave functions can be obtained from numerical calculations of the four-quark form factors and quasi TMDPDFs and TMDWFs on the lattice [26,27]. On the other hand, one can also make use of the QCD factorization to obtain the CS kernel from quasi TMDPDFs and TMDWFs. The first results for the CS kernel based on these proposals have been published recently [29][30][31][32][33]. The quasi TMD-WFs approach to the CS kernel requires two-point function calculations and potentially can reach the light-cone limit with relatively small hadron momenta.
In this work, we present a state-of-the-art calculation of the CS kernel, based on a lattice QCD analysis of quasi-TMDWFs with N f = 2 + 1 + 1 valence clover fermions on a staggered quark sea with one-loop matching accuracy. A single ensemble with lattice spacing a ≃ 0.12 fm, volume n 3 s ×n t = 48 3 ×64, and physical sea-quark masses is used. In order to improve the signal-to-noise ratio, we tune the light-valence quark masses such that m π = 670 MeV. The CS kernel is then extracted through the ratios of the quasi-TMDWFs and the perturbative matching kernels at different momenta, P z = 2π/n s × {8, 10, 12} = {1.72, 2.15, 2.58}GeV. This corresponds to Lorentz boost factors γ = {2.57, 3.21, 3.85}, respectively. This analysis improves the previous ones [30,32] by taking into account the one-loop perturbative contributions, and by analyzing systematic uncertainties from operator mixing, higher-order corrections from the scale dependence, and higher power corrections in terms of 1/P z .
The remainder of this paper is organized as follows. In Sec. II, we present the theoretical framework to extract the CS kernel from quasi-TMDWFs. Numerical results for quasi-TMDWFs and CS kernel are presented in Sec. III. A brief summary of this work is given in Sec. IV. More details about the analysis are collected in the appendix.

II. THEORETICAL FRAMEWORK
In this section, we review the necessary theoretical background for the present calculation. We present the definitions of CS kernel and rapidity evolution, and introduce the quasi-TMD wave functions. We then discuss the factorization of the quasi-TMDWFs and its connection with the CS kernel.

A. Collins-Soper Kernel and Rapidity Evolution
Unlike the collinear lightcone PDFs and distribution amplitudes, the TMDPDFs and TMDWFs depend on both the renormalization scale µ and an additional rapidity renormalization scale. The latter arises because the matrix elements also suffer from so-called rapidity divergences that require a dedicated regulator [1,34,35]. In TMD factorizations, the contributions of hard, i.e. highly offshell, modes to the tree process are usually calculated in the dimensional regularization scheme. Collinear modes, which are related to highly-boosted partons in distinct directions, and soft modes, whose typical momentum are at the order Λ QCD , share the same virtuality, and are only distinguishable by their rapidity. In calculations using regularization schemes such as dimensional regularization, which only regulate ultra-violet divergences, one will encounter additional rapidity divergences that arise in soft and collinear matrix elements when integrating over rapidity, and have to be resolved using a dedicated regulator. After the later regularization, TMDPDFs and TMDWFs acquire an additional rapidity scale dependence. This dependence should cancel in theoretical predictions for physical observables.
The CS kernel K (b ⊥ , µ), known as the rapidity anomalous dimension, encodes the rapidity dependence of the TMD distributions [1,2]: where f TMD denotes any leading twist TMDPDF or TMDWF. The TMD distributions depend on the longitu-dinal momentum fraction x, transverse spatial separation b ⊥ , which is the Fourier-conjugate to the transverse momentum k ⊥ , as well as the renormalization scale µ and rapidity scale ζ which is related to the hadron momentum. The µ-dependence of CS kernel K (b ⊥ , µ) satisfies the renormalization group equation (RGE): Here Γ cusp (α s ) = α s C F /π+O(α 2 s ) is the cusp anomalous dimension, which has been calculated in perturbation theory up to 2-loop in Ref. [36], and 3-loop in Ref. [37]. The solution to the RGE can be expressed as: For large b ⊥ with b −1 ⊥ Λ QCD , the CS kernel becomes nonperturbative, which is represented by the non-cusp anomalous dimension K (α s (1/b ⊥ )) in the above equation.
In the past decades, the CS kernel has been widely studied in global fits of TMD parton distributions [12][13][14][15][16][17][18][19][20]. The explicit form in the nonperturbative region can only be parametrized by extending the perturbative expressions at small b ⊥ , which inevitably introduces systematic uncertainties. A direct calculation of TMD-PDFs and the relevant CS kernel on the lattice was an almost insurmountable hurdle until the establishment of LaMET [21,22]. A remarkable recent development in LaMET is that these quantities can be accessed through the corresponding quasi observables [25][26][27][28].

B. Quasi TMD Wave Functions
As stated above, one can define the quasi TMDWFs for a highly-boosted pseudoscalar meson along the zdirection with large momentum P z as: where x r = x− 1 2 . The unsubtracted quasi TMDWFΦ ±0 is defined as an equal-time correlator containing nonlocal quark bilinear operator with staple-shaped gauge link: For a pseudoscalar mesonic state, the Dirac structure Γ can be chosen as γ z γ 5 or γ t γ 5 , which approaches the leading-twist structure γ + γ 5 in the light-cone limit. With a large but finite P z , the difference between γ z γ 5 and γ t γ 5 is suppressed by powers of 1/P z . Technically, one can also use a combination of them, such as (γ z γ 5 + γ t γ 5 ) /2 to minimize power corrections, and more details can be found in Sec. III D. Various combinations were also explored in Ref. [32]. The superscript "0" inΦ ±0 indicates bare quantities. The linear divergences come from the self-energy of the gauge-link, and does not appear as a pole at d = 4 in dimensional regularization. The Euclidean gauge link in U ⊐,±L is defined as where ξ z = −ξ · n z . The ±L corresponds to the farthest position that the gauge link can reach in positive or negative n z direction on a finite Euclidean lattice. This is depicted as the blue and red lines in Fig. 1.
Since the linear divergence is associated with the gauge-link, it can be removed by a similar gauge-link with the same total length. An optional choice is to make use the Wilson loop, denoted as Z E . The Wilson loop can be chosen as the vacuum expectation of a flat rectangular Euclidean Wilson-loop in the z-⊥ plane: Here the length of Z E is twice that of the staple-shaped gauge-link U ( †) z in the z direction, and thus it is anticipated that the square root of Z E (2L, b ⊥ , µ) cancels the linear divergence and heavy quark potential in the gaugelink. There are residual logarithmic divergences from the vertex of Wilson line and light quark, which can be renormalized in dimensional regularization [38]. As these logarithmic divergences are independent of z, b ⊥ , P z and L, they will explicitly cancel out when the ratio of quasi TMDWFs is studied.

C. Factorization of Quasi TMDWFs
With the help of soft function, the infrared contributions in the subtracted quasi TMDWFs can be properly accounted for such that the infrared structures for the quasi TMDWFs and light-cone ones are matched. This implies a multiplicative factorization theorem in the framework of LaMET [25][26][27][28]39]: where the superscript ± in Eq. (9) corresponds to the direction in the Wilson line, Ψ ± is the TMDWFs defined in the infinite momentum frame. The reduced soft function S 1/2 r (b ⊥ , µ) emerges from the different soft gluon radiation effects inΨ ± and Ψ ± [28]. The mismatch of the rapidity scale ζ and ζ z can be compensated by the CS kernel K(b ⊥ , µ). Both S and K are independent of the ± choice. H ± is the 1-loop perturbative matching kernel [28]: With the abbreviations ℓ ± = ln (−ζ z ± iǫ)/µ 2 , and It should be noticed that H ± contains nonzero imaginary parts in ℓ ± and ℓ ± . While the imaginary parts in ℓ ± /l ± are constants, the ones in the double logarithms ℓ 2 ± /l 2 ± are momentumdependent.
A characteristic behavior of Eq.(9) is that this factorization is multiplicative [24], which indicates that hard gluon contributions are local. This is due to the fact that hard gluon exchange between the quark and antiquark sectors in quasi TMDWFs is power suppressed: if there were such a hard gluon, the spatial separation between its attachments is much smaller than b ⊥ , resulting in power suppression compared to the typical hard Leading-power reduced graph for pseudoscalar meson quasi TMDWFs. Here C, S denote the collinear and soft sectors of the infra-red structure, while H denote the hard contributions. Since the hard-gluon exchange between the quark and anti-quark is power suppressed, the hard parts are disconnected with each others, and therefore the factorization of quasi-WF amplitude is multiplicative. mode contributions. Thus, at leading power the factorization of quasi TMDWFs are multiplicative. This feature is illustrated in Fig.2, in which the collinear, soft and hard sub-diagrams represent the pertinent contributions. Further, ζ z arising from Lorentz-invariant combinations of collinear momentum modes, will provide the natural hard scale of the hard sub-diagram. More detailed explanations for the factorization of quasi TMDPDFs in LaMET are given in the recent review [27].

D. Collins-Soper Kernel From Quasi TMDWFs
From Eq.(9), one can see that the momentum dependence in quasi TMDWFs provides an option to determine the CS kernel. This can be written in a way similar to Eq.(1) [28]: where K (b ⊥ , µ) denotes the same kernel as in Eq.(1), and does not depend on the hard scale ζ z for large P z . Unlike TMDWFs, the quasi distributions also contain hard contributions, whose rapidity dependence is represented by the perturbative G ± as a function of hard scale ζ z . From the ζ z dependence of quasi TMDWFs, we can see that when P z → ∞, the large logarithms in P z are partially absorbed into K (b ⊥ , µ) and the remanent is incorporated in the perturbative matching kernel. Therefore, both the matching kernel H ± and an exponential of CS kernel K (b ⊥ , µ) are needed to describe the dependence on ζ z of quasi TMDWFs. In order to extract the CS kernel K (b ⊥ , µ) explicitly, one can make use of Eq.(9) with two different large momenta P z 1 = P z 2 ≫ 1/b ⊥ but the same scale ζ z . Taking a ratio of these two quantities gives where the reduced soft function S r (b ⊥ , µ) and TMDWFs Ψ ± (x, b ⊥ , µ, ζ) have been canceled in the ratio. Therefore, the CS kernel K (b ⊥ , µ) can be extracted through .
Note that the extracted result is formally independent of x and P z 1/2 at leading power, and bothΨ + andΨ − can be used to extract K (b ⊥ , µ). This is derived at the leading power in the factorization scheme and might be undermined by power corrections. Accordingly, in order to reduce the systematic uncertainties, we take the average: .
The details will be discussed in Sec. III E.

III. NUMERICAL SIMULATIONS AND RESULTS
In this section, we present our lattice QCD results. We start with the lattice setup, followed by results for quasi TMDWFs with two-point correlations. The Wilson loop results are discussed in subsection III.C. Subsection III.D studies the operator mixing effects. Our main result on CS kernel is presented in subsection E. The final subsection E includes some overall discussions.

A. Lattice setup
Our numerical simulations use N f = 2 + 1 + 1 valence clover fermions on a highly improved staggered quark (HISQ) sea [40] and a 1-loop Symanzik improved gauge action [41], generated by the MILC collaboration [42] using periodic boundary conditions. In the calculations, we use a single ensemble with the lattice spacing a ≃ 0.12 fm and the volume n 3 s × n t = 48 3 × 64 at physical sea-quark masses. In order to increase the signal-to-noise ratio, we tune the light-valence quark masses to the strange-quark one, namely m sea π = 130 MeV and m val π = 670 MeV, which could generate some non-unitarity effects. On the other hand, the Collins-Soper kernel only depends weakly on quark mass, and we may consider valence quarks are strange-like, namely the hadrons involved are kaons.

B. Quasi TMDWFs From Two-Point Correlators
In order to calculate the quasi TMDWFs defined in Eq.(4), we generate Coulomb-gauge wall-source propagators, where (t ′ , y) and (t, x) denote the space-time positions of source and sink. Then one can construct the two-point function (2pt) related to the quasi TMDWFs in Eq. (4): The quark momentum p = ( 0 ⊥ , p z ) is along the z-direction, and each of two quarks carries half of the hadron momentum. Thereby the hadron momentum satisfies P = 2 p. The anti-quark propagator can be obtained from Eq.(10) by applying γ 5 -hermiticity S w (x, y) = γ 5 S w (y, x) † γ 5 . As mentioned above, the Dirac structures are chosen as Γ = γ z γ 5 and γ t γ 5 that can be projected onto the leading twist light-cone contributions in the large P z limit. By generating the wall-source propagators with quark momenta p z = ±{4, 5, 6} × 2π/(n s a), and three segments of gauge links following Eq.(6), one can construct the two-point correlation functions on the lattice. With the help of reduction formulas, the C ± 2 can be parametrized as where A w (p z ) is the matrix element of the pseudoscalar meson interpolating field with Coulomb gauge fixed wall source and A p is the one for a point source (sink). These terms as well as the factor E −1 = 1/ m 2 π + (P z ) 2 are cancelled by the local two-point function C ± 2 (0, 0, P z ; p z , 0, t) at the same time slice. Thus the remaining ground-state matrix elementΦ ±0 is normalized. The ratio of nonlocal and local two-point Comparison of two-state fit and one-state fit to extractΦ ±0 (z, b ⊥ , P z , L).
Taking {z, b ⊥ , P z , L} = {0a, 2a, 24π/ns, 6a} as example, we can see that the twostate fit works for t ∈ [2a, 8a] while the one-state fit works for t ∈ [5a, 8a]. The fitted results are consistent with each other, while the one-state fit is more conservative.
functions can be parametrized as where C 2 in the denominator is a local correlator.
In the above parametrization, the excited-state contributions are collected into the c 0 term, and ∆E denotes the mass gap between the ground and first excited state. With the increase of Euclidean time, contributions from the excited state decay and the plateau obtained for R ± (z, b ⊥ , P z , L, t) at large times reflects the ground-state contributionΦ ±0 . We employ two methods to extractΦ ±0 , namely the two-state fit directly using Eq. (18), and the one-state fit by setting c 0 = 0. With a large enough Euclidean time, Fig. 3 exhibits a comparison using two methods for the case with small {z, b ⊥ }. From this figure, one can see that the one-state fit result is consistent with two-state fit one but gives a more conservative error estimate. The two-state fit works at t ∈ [2a, 8a] and the one-state fit works at the plateau region t ∈ [5a, 8a]. While the excited state contamination would dominate for the two-state fit in a very high precision especially for the cases with large {z, b ⊥ }. So with current accuracy, we adopt the more conservative results from the one-state fit in the following analysis. More details can be found in Appendix A.

C. Wilson Loop Renormalization
The unsubtracted quasi TMDWF matrix elements Φ ±0 (z, b ⊥ , P z , a, L) extracted from the joint fit of the two-point function contain a factor e −δm(2L+b ⊥ ) from the linear divergence, the heavy quark effective potential factor e −V (b ⊥ )L and logarithmic divergences Z O : where Z O has logarithmic dependence on lattice spacing a.
The linear divergence in e −δm(2L+b ⊥ ) comes from the self-energy of the Wilson line [44][45][46][47], where δm contains a term proportional to 1/a and a non-perturbative renormalon contribution m 0 : Note that the exponent of the linear divergence term is proportional to the total length of the Wilson link, e.g. 2L + b ⊥ for the staple link. Due to this factor, the numerical value for a Wilson loop dramatically decreases for small a and large L.
The heavy quark effective potential term e −V (b ⊥ )L comes from interactions between the two Wilson lines along the z direction in the staple link. The heavy quark effective potential V (b ⊥ ) is often used to determine the lattice spacing of an ensemble.
The logarithmic divergence Z O comes from the vertices involving the Wilson line and light quark. The logarithmic divergence up to leading order, resumed by renormalization group equation, and matched to MS scheme is [48,49] where Λ latt QCD is different from that in MS scheme. One can use both to effectively absorb higher order contributions [50].
In this work, the Wilson loop renormalization method [46,[51][52][53][54] is adopted, in which the Wilson loop Z E defined in Eq.(8) contains linear divergence and heavy quark effective potential: According to Ref. [49], δm in the Wilson loop is the same as that in hadron matrix elements, and thus it is anticipated that the linear divergence is removed when dividing by Z E (2L, b ⊥ , a): As shown in Fig. 4, the subtracted quasi TMDWFs tend to be a constant when L ≥ 0.4 fm. We then use the subtracted quasi TMDWFs defined as However, it is anticipated that there is a residual logarithmic divergence Z O : In the extraction of the CS kernel, a ratio of quasi TMD-WFs is adopted and accordingly the residual logarithmic divergence Z O is canceled.

D. Higher-Twist Effects in Operators
For a pseudoscalar meson on a Euclidean lattice, both Γ = γ t γ 5 and γ z γ 5 project onto the leading twist lightcone distribution amplitude, i.e. γ + γ 5 in the arge P z limit. The differences between them arises from power corrections in terms of M 2 / (P z ) 2 .   5 shows the comparison of the λ = zP zdependence of quasi TMDWFs with Γ = γ t γ 5 and γ z γ 5 at P z = 24π/n s ≃ 2.58 GeV. It can be seen from the plots that there are some differences between the two sets of results for the real part in the small λ region. The differences are expected to decrease with increasing P z , and the correlators with γ t γ 5 and γ z γ 5 will gradually converge to the light-cone from opposite directions. Besides, in light-cone coordinate, γ t and γ z can be represented by γ + and γ − , Due to the momentum along the light-cone, operators with γ − γ 5 correspond to higher order terms of TMDWFs. Therefore, power corrections arising from finite P z are likely to be eliminated in the average of these two terms: For a quantitative analysis see the appendix of [30]. The operator mixing effect reaches order 5% [30], which is much smaller than the systematic uncertainties discussed in next subsection. According to our numerical simulations, the subtracted quasi TMDWFs in coordinate spaceΦ ± (z, b ⊥ , P z ) as a function of λ = zP z are complex, which is shown in Fig. 6. The examples are the real and imaginary part ofΨ ± (x, b ⊥ , P z ) with P z = 24π/n s , b ⊥ = 2a and 4a. To determine quasi TMDWFs in momentum spacẽ Ψ ± (x, b ⊥ , P z ), we use a Fourier transformation (FT) Due to the imaginary part ofΦ ± (z, b ⊥ , P z ), Ψ ± (x, b ⊥ , P z ) also has an imaginary part. We obtain the quasi TMDWFs in momentum space for both real and imaginary part ofΨ ± (x, b ⊥ , P z ) shown in Fig.7, by taking P z = 24π/n s and b ⊥ = 2a and 4a as examples. We truncate the FT at z min and z max . The deviation ofΦ ± (z min , b ⊥ , P z ) andΦ ± (z max , b ⊥ , P z ) from zero is a measure of the resulting truncation error. For the largest range of z values we could realize numerically, z min = −1.44fm, z max = 1.44fm, this error is still noticeable. This brute-force truncation of the FT leads to an oscillatory behavior of TMDWFs. This oscillation inΦ(x, b ⊥ , P z ) can be eliminated by an appropriate extrapolation forΦ(z, b ⊥ , P z ) as a function of zP z before Fourier transformation. While the signal-to-noise ratio of our data is not smooth enough, the brute-force Fourier transformation is adopted.

E. Collins-Soper Kernel From Quasi TMDWFs
The CS kernel governs the rapidity evolution and thus is independent of the momentum fraction of the involved parton. But as indicated in Eq.(9), the factorization formula works only when xP z ≫ Λ QCD , and could be invalid in the end-point regions x → 0, 1. Power corrections are presumably of the form 1/ (xP z ) 2 or 1/ (xP z ) 2 . Therefore, the numerical CS kernel is fitted by a function of x, P z 1 and P z 2 and is written as K(b ⊥ , µ, x, P z 1 , P z 2 ), Here K(b ⊥ , µ, x, P z 1 , P z 2 ) are extracted from the perturbative matching kernels and quasi TMDWFs using 1-loop matching. They will have power corrections of teh form O 1/ (xP z ) 2 and O 1/ (xP z ) 2 . In order to extract the leading power contributions, we adopt the following parametrization where A is the coefficient accounting for the leading higher power contributions, and can be determined through a joint fit of different lattice data in the regions not so close to x = 0, 1. Fig. 8 presents the physical CS kernel K(b ⊥ , µ) with b ⊥ = {0.12, 0.24, 0.36, 0.48, 0.60}fm. By employing three cases of quasi TMDWFs with P z = {8, 10, 12} × 2π/n s , one can extract K(b ⊥ , µ, x, P z 1 , P z 2 ) with P z 1 /P z 2 = 10/8 and 12/8, shown as the different colored bands. Except in the end-point regions (x < 0.2 or x > 0.8), the lattice data is flat and reflects the leading power contribution, which conforms with expectations. Using the parametrization formula Eq.(30), the physical CS kernel K(b ⊥ , µ) can be determined by fitting the data, shown as the green band. As mentioned before, at large b ⊥ , the quasi TMDWFs show oscillations due to the truncation of the Fourier transformation, which also affect the extracted K(b ⊥ , µ, x, P z 1 , P z 2 ), as shown in the lower panel of Fig.8. This oscillation effect can be in principle removed once larger z data becomes possible, or if one knows how to extrapolated the current data to the larger z or λ region.
Theoretically the physical CS kernel is purely real, however, there still exists a residual imaginary part at 1-loop matching. As discussed above, this imaginary term comes from the perturbative matching kernel. It is easily to prove that H ± (zP z 2 , µ)/H ± (zP z 1 , µ) contains imaginary part while the lattice results for Ψ ± (x, b ⊥ , µ, P z 1 )/Ψ ± (x, b ⊥ , µ, P z 2 ) are nearly real, shown as 9. Therefore, we consider this imaginary part as systematic uncertainty from factorization theorem. This uncertainty can be expressed as where Im [K + (b ⊥ , µ)] represents the numerical imaginary part of extracting K(b ⊥ , µ) only byΨ + : .
It should be noticed that the perturbative matching kernel H + (zP z , µ) is the complex conjugate of H − (zP z , µ), that is the imaginary parts in these terms can be cancelled each other when we employ the average of H + (zP z 2 , µ)/H + (zP z 1 , µ) and H − (zP z 2 , µ)/H − (zP z 1 , µ). Therefore, as the final result, ] /2 to reserve the real part, and regard the imaginary contributions as our systematic uncertainty.

F. Results and Discussions
One should notice that the Wilson loop renormalized quasi TMDWF on the lattice (Eq.(24)) has a scale dependence on a. If one converts it to the MS scheme through dividing it by Z O (Eq.(21)), the scale µ is introduced. In principle, one should convert the Wilson loop renormal- ized quasi TMDWF to the MS scheme since our factorization formula works there. However, since Z O has no dependence on momentum P z , it cancels in the ratio of quasi TMDWFs, so does the scale dependence. So, one does not need to do the scheme and scale conversion of the quasi TMDWF during the extraction of CS kernel.
The extracted CS kernel from the combined fit of the ratios of quasi TMDWFs with different momenta are shown as the red data points in Fig. 10. In this figure we exhibit two kinds of errors for K(b ⊥ , µ), in which the smaller ones denote statistical uncertainties while the larger ones include both statistical and systematical uncertainties. In the small-b ⊥ region, systematical uncertainties are dominant due to the large power and the nonzero imaginary part.
As a comparison, we also give the tree-level matching result for the CS kernel. With the leading order matching kernel H(xP z , µ) = 1 + O(α s ), Eq. (13) simplifies to the ratio of quasi TMDWFsΦ at z = 0 with momentum P z 1 /P z 2 . The blue dots in Fig. 10 denote the results obtained for tree-level matching, for which only statistical uncertainties are shown.
We compare our results with the ones from perturbative calculations, phenomenological extractions as well as the lattice results determined by other collaborations.
The black solid and dashed lines in the upper panel of Fig. 10 indicate the perturbative results up to 3-loops, with a running coupling constant α s (µ = 1/b ⊥ ). The perturbative calculations work well in small b ⊥ region (b ⊥ ≪ 1/Λ QCD ), while will diverge with b ⊥ increasing. In contrast, the lattice calculation will give accurate predictions in the nonperturvative region, while due to the power corrections, it might suffer large systematic uncertainties in small b ⊥ region.
Similar to this work, the results of LPC [30] and ETMC/PKU [32] are also extracted from quasi TMD- , with b ⊥ = 3a as a function of momentum fraction x. The imaginary parts of both cases are close to zero. WFs through a tree-level matching. Adopting the oneloop matching formula, as well as considering the operator mixing effects will help ones to reduce the systematic uncertainties and obtain more precise results. In addition, considering the different directions of gauge link will help us to eliminate the contributions from unphysical imaginary part, and then improve the accuracy of our results.
In another way, the SWZ [33] and SVZES [31] results are obtained from quasi-TMDPDFs. Compared with the complicated nucleon correlation functions, the meson ones are much easier to obtain better signals. Besides, the wave functions of meson are nearly symmetric in x-space, it thereby more convenient to parametrize the oscillation effects and obtain the physical results in large P z limit like Fig.(8). In addition, the light meson is more easily to reach a larger boosted factor, one can see from the small b ⊥ region, the results from quasi TMDWFs fit well with the perturbative calculations than the ones from quasi TMDPDFs.
The lower panel of Fig.(10) shows the comparison with phenomenological results. SV19 [19] and SIYY15 [15] use a parameterization with perturbative and nonperturbative parts. However Pavia19 [20] obtained their result with the factorization of TMDPDFs, obtaining the CS kernel from the rapidity derivative. In addition they fit parameters from the Drell-Yan data to obtain their phenomenological CS kernel. The results from different methods exhibit obviously inconsistency in the nonperturbative region. Our result shows a better consistency with SV19. FIG. 10. The upper panel shows the comparison of our results K(b ⊥ , µ) and K0(b ⊥ , µ) with the lattice calculations by SWZ [33], LPC [30], ETMC/PKU [32] and SVZES [31], as well as the perturbative calculations up to 3-loop. K(b ⊥ , µ) denotes the CS kernel extracted through 1-loop matching, and whose uncertainties correspond to the statistical errors and the systematic ones from the non-zero imaginary part. K0(b ⊥ , µ) denotes our tree-level results, only with statistical uncertainties. The lower panel shows the comparison of our result with phenomenological extractions: SV19 [19], Pavia19 [20] and SIYY15 [15] give phenomenological parameterizations of CS kernel fitted to data from high energy collision processes like Drell-Yan.

IV. SUMMARY AND OUTLOOK
In this work, we have calculated the CS kernel on a MILC lattice configuration in the large momentum effective theory framework. Comparing with our previous studied [30], the one-loop matching kernel has been adopted in this study, and several hadron momenta were used to extract the CS kernel. We found that in the small b ⊥ region, our results are consistent with perturbative QCD. In large b ⊥ region, our results seem consistent with other lattice calculations in the literature within uncertainties.
For our future studies, we need to use lattice configurations with multiple lattice spacing to understand the finite lattice spacing effects. We would use a valence quark mass consistent with the sea quark one to reduce the nonunitarity effects. One such effect might be the imaginary part of the meson wave function which seems inconsistent with perturbative calculation at present time. Clearly, all these explorations will take more computational resources.