Transverse-momentum-dependent wave functions and Soft functions at one-loop in Large Momentum Effective Theory

In large-momentum effective theory (LaMET), the transverse-momentum-dependent (TMD) light-front wave functions and soft functions can be extracted from the simulation of a four-quark form factor and equal-time correlation functions. In this work, using expansion by regions we provide a one-loop proof of TMD factorization of the form factor. For the one-loop validation, we also present a detailed calculation of ${\cal O}(\alpha_s)$ perturbative corrections to these quantities, in which we adopt a modern technique for the calculation of TMD form factor based the integration by part and differential equation. The one-loop hard functions are then extracted. Using lattice data from Lattice Parton Collaboration on quasi-TMDWFs, we estimate the effects from the one-loop matching kernel and find that the perturbative corrections depend on the operator to define the form factor, but are less sensitive to the transverse separation. These results will be helpful to precisely extract the soft functions and TMD wave functions from the first-principle in future.

In large-momentum effective theory (LaMET), the transverse-momentum-dependent (TMD) light-front wave functions and soft functions can be extracted from the simulation of a four-quark form factor and equal-time correlation functions. In this work, using expansion by regions we provide a one-loop proof of TMD factorization of the form factor. For the one-loop validation, we also present a detailed calculation of O(αs) perturbative corrections to these quantities, in which we adopt a modern technique for the calculation of TMD form factor based the integration by part and differential equation. The one-loop hard functions are then extracted. Using lattice data from Lattice Parton Collaboration on quasi-TMDWFs, we estimate the effects from the one-loop matching kernel and find that the perturbative corrections depend on the operator to define the form factor, but are less sensitive to the transverse separation. These results will be helpful to precisely extract the soft functions and TMD wave functions from the first-principle in future.

I. INTRODUCTION
The exploration of underlying structures of hadrons has always been one of the most important frontiers in particle and nuclear physics. The one-dimensional and three-dimensional light-front wave functions (LFWFs) are important physical quantities describing the distributions of constituents' momentum in the hadron, and reflect the non-perturbative internal structure of hadrons [1][2][3]. As for a light Nambu-Goldstone boson, the LFWFs also help us to understand the chiral symmetry breaking [4][5][6][7].
It was firstly noticed that in an exclusive process [1,2] the non-perturbative LFWFs for a given Fock state are required. As an inevitable input, LFWFs also play an important role in theoretical analyses of B meson weak decays [8][9][10][11][12], which are of great values for the test of the standard model (SM) and the search for new physics beyond the SM. The unprecedented high precision of measurements at the current and forthcoming experimental facilities strongly request the improvement of the accuracy of these non-perturbative physical quantities from quantum chromodynamics (QCD).
In 1970s, the formation of one-dimensional hadronic wave function, namely, light-cone distribution amplitudes (LCDAs) was established in the large momentum limit under light-front quantization [1]. Results for a few lowest moments of the LCDAs were firstly obtained from QCD sum rules [13,14], and since then many progresses have been made in extracting the moments of LCDAs in the past decades. Despite of these progresses, a complete knowledge of meson wave functions from the first principle is not well-established yet.
An obvious difficulty in calculating LCDAs lies in the * Corresponding author: zengj@sjtu.edu.cn fact that it is inherently non-perturbative and needs to be treated by methods such as lattice QCD (LQCD). However, LCDAs belong to light-cone correlations of quark/gluon field operators, and thus contain explicit time dependence, which cannot be directly calculated by lattice field theory defined in the Euclidean spacetime. In this respect, only the moments of LCDAs, namely matrix elements of the local operators, can be performed in the traditional LQCD approach [15][16][17][18][19][20][21].
A remarkable approach to circumvent the above problem is proposed in Ref. [22,23], which is now systematically formulated as large momentum effective theory (LaMET). In LaMET, one can construct the directly computable hadron matrix elements with non-local operators, named as quasi-distributions, on the lattice. Through a perturbative matching, the corresponding LCDAs can be accessed [22][23][24][25]. Many inspiring results on LCDAs were reported in recent years, and reviews of recent developments can be found in Refs. [24,25].
Compared with LCDAs and parton distribution functions (PDFs), transverse-momentum-dependent wave functions (TMDWFs) and TMD parton distribution functions (TMDPDFs) provide more versatile information on the internal three-dimensional structure of hadrons, which are also relevant for observables with transverse momentum dependent (TMD) distributions of final-state particles in high-energy experiments. For instance, TMDWFs have been applied to calculate various transition form factors such as the pion electromagnetic form factor [26,27], the proton form factors [28][29][30][31] and exclusive B decays [8,32]. Therefore, it is highly prerequisite to further investigate the three-dimensional LFWFs from the first-principle QCD.
A very important progress in LaMET is that TMD distributions can be accessible through the Euclidean equal-time correlations [33][34][35][36][37]. In Ref. [35], it has been demonstrated that the form factor of a bi-local fourquark operator, calculable on the lattice, can be fac-torized into TMDWFs, a universal soft factor (function) and the matching kernel through QCD factorization at large momentum transfer. A combined analysis of quasi-TMDWFs on lattice allows a direct extraction of the universal soft function and TMDWFs. Based on these proposals, Lattice determinations of the rapidity evolution anomalous dimension, namely Collins-Soper (CS) kernel [38], can be found in Refs. [39][40][41][42][43][44].
It is anticipated that in the large momentum limit, the form factor can be expressed as a convolution of TMD-WFs, soft functions and a hard kernel. In this work, we aim to present a complete one-loop analysis of these quantities and provide the necessary details to the proposal in Refs. [35,37]. We use the expansion by regions and provide a proof of the TMD factorization for the form factor. In the explicit calculation factor, we adopt a modern technique based on the integration by part (IBP). With these results, we will demonstrate the cancellation of the infrared divergences, and explicitly validate the TMD factorization. Finally, we extract the hard kernels for the form factor and quasi-TMDWFs through the factorization, which will be useful for a precision determination of TMDWFs and soft functions. We also make use of the lattice data from Lattice Parton Collaboration on quasi-TMDWFs [44] and show that the perturbative corrections to soft functions depend on the Lorentz structures, and the magnitude can reach the order (10 − 30)%. The results are found to be less sensitive to the transverse separation, which may imply a factorized form for the quasi-TMDWFs. As a comparison, we also give a phenomenological parametrization which contains explicit dependence on the transverse separation.
The rest of this work is organized as follows. In Sec. II we will briefly introduce the concept of TMDWFs. In Sec. III, the one-loop perturbative results for TMDWFs, soft functions, form factors and Wilson loops will be presented in order. To regularize the rapidity divergence in TMDWFs and soft function, the delta regulator will be used in the calculation. Based on these results, we will explicitly validate the TMD factorization of four quark form factors and extract the short-distance hard kernel at one-loop level. In Sec. IV, the one-loop perturbative results for quasi-TMDWFs will be presented. In Sec. V, we use the hard kernel and calculate the effects to extract the soft functions. We conclude this work in Sec. VI. In the appendix, we collect some details in the calculation.

II. TMD WAVE FUNCTIONS
In this section, we will follow the spirit of Ref. [37] and give a self-contained description of TMDWFs.
An intuitive identification of LFWFs is the light front correlation functions between hadron state and QCD vacuum, in which light-like gauge links extending to infinities are required to maintain the gauge invariance. This allows the identification of LF divergences as rapidity divergences, known in the literature of TMD physics.
In high energy limit, light quarks and gluons inside a hadron move on the lightcone, and are generally named as partons. In parton physics, light-front quantization (LFQ) is a useful formalism to handle the hadron states. It provides a Hamiltonian description of QCD similar with the diagonalized Hamiltonian in non-relativistic quantum mechanics aŝ where |Ψ n denotes a QCD bound state [3], and P + = (P 0 + P z )/ √ 2. The wave functions obtained in this picture can in principle be used to calculate all partonic densities and correlations functions.
In the infinity momentum frame (IMF), making an IR cut-off on the longitudinal momentum scale k + = ǫ and taking all physics below it into renormalization constants, one can get an effective Hilbert space and obtain an effective LF theory with trivial vacuum, Here |0 is the vacuum of LFQ, and a i denotes the annihilation operator of all kinds of possible partons. Therefore, in the LF gauge A + = 0 the hadron can be expanded in terms of the superposition of all kinds of possible Fock states [3], where a † i is the creation operator of partons on the lightfront, and the phase-space integral takes a light-core decomposition dΓ n = dk + d 2 k ⊥ 2k + (2π) 3 . The ψ n (x i , k i⊥ ) is the LFWF, where x i denotes the set of momentum fractions of each parton, and k i⊥ is the corresponding transverse momentum. The summation over index n sums all possible partons of hadron state, and the multiplication over index i multiplies all partons which that in the state.
With the truncation k + ≥ ǫ, we can write the above expansion in form of invariant matrix elements, With the inclusion of gauge-invariance and regularizations, this invariant matrix element will become the correlator matrix elements.
For a generic notation φ i denotes all kinds of partons including the quark fields ψ and gluon fields A µ , with the index 'i' to label the field. We introduce a gaugeinvariant field Φ i which contains gauge-link along the light-cone direction n, pointing to positive or negative infinity: with the light-like Wilson line W ± n (ξ) where P is a path order. Then the generic naive threedimensional LFWFs, namely TMDWFs, are written as In the above N k=0 x k = 1, N k=0 λ k = 0, and each x i≥1 are longitudinal momentum fractions carried by partons satisfying 0 < x k < 1. Likewise, when P has no transverse component, the transverse coordinate b ⊥ can be shifted by an overall constant without any effect.
The ultraviolet and infrared divergences from these amplitudes can be regulated in dimensional regularization (DR) with modified minimal subtraction (MS)scheme. However, in TMDWFs, there is a new type of divergence from light-like gauge-links extending to infinities, which is called rapidity divergence. It arises from the collinear gluons radiation with the momentum fraction approaching zero but cannot be regulated by dimensional regularization. There are multiple methods to regulate the rapidity divergence, and an option is the so-called δ regulator [45,46]. In this scheme, the gauge-link is modified as [37] where δ − is a positive quantity to characterize the rapidity divergence. It breaks the gauge-invariance, but the breaking effects approach zero when δ − → 0. The regularization for the other light-cone direction is similar with the regulator δ + . Then the TMDWFs are written as where the fields Φ i are now defined as As δ − → 0, TMDWFs and TMDPDFs diverge logarithmically, and the remanent finite part also depends on the rapidity regulator. Therefore, the naive TMDWFs in Eq. (9) can not solely absorb all nonperturbative dynamics in the factorization for physical observables. One must remove all divergences and rapidity regularization scheme dependencies in ψ, in a way similar with removing UV divergences in physical quantities. These are accomplished with the help of soft functions to be introduced in the next section.
For a pseudoscalar pion, the TMDWF is defined as where b µ = (0, b ⊥ , 0) is the transverse space coordinate of the light-quark field, f π is the decay constant of pion. The factor 1 −ifπ P + comes from a normalization in terms of a hadronic local operator matrix element,

III. TMDWFS IN LAMET
At present, a most systematic approach to solve nonperturbative QCD is lattice field theory [47]. Though quantities on the lightcone can not be straightforwardly implemented on the lattice, LaMET offers a practical way to carry out the program of light-front quantization (LFQ) [37]. In a certain sense, the quantization using tilted light-cone coordinates [48] is similar to the spirit of LaMET [37]. Therefore, a practical implementation of LaMET can be done through lattice calculations. While LFQ may provide an attractive physical picture for the proton, the Euclidean equal-time formulation is more practical for carrying out the calculations, and LaMET serves to bridge them.
The relation between partonic observables on the LF and the properties of a hadron with a large momentum is not one to one. There are infinite possible Euclidean operators in the large-momentum proton that generate the same LF observable. This is because the large-momentum physical states have built-in collinear (as well as soft) parton modes, and upon acting on a Euclidean operator they help to project out the leading LF physics. All operators projecting out the same LF physics form a universality class. In the operator formulation for parton physics such as soft collinear effctive theory (SCET) [49], one uses LF operators to project out parton physics off the external states of any momentum, including P = 0. Concepts such as the universality class have been explored in critical phenomena in condensed matter physics, where systems with different microscopic Hamiltonians can have the same scaling properties near their critical points. Critical phenomena correspond to the infrared fixed points of the scale transformation and are dominated by physics at long-distance scales. In this case, parton physics arises from the infinite momentum limit, which is a ultraviolet fixed point of the momentum renormalization equations (RGEs). It is the longitudinal short-distance physics that is relevant at the fixed point. However, the short distance here does not mean that everything is perturbative. The part that is nonperturbative characterizes the partonic structure of the meson. The critical region P → ∞ acts as a filter to select only the physics that is relevant, so universality classes emerge. It has been pointed out that the soft function can be obtained from a form factor of a pseudoscalar light-meson state [35]: where ψ a,b,c,d are light quark fields of different flavors.
The factor 1 f 2 π P1·P2 comes from the normalization of two local hadronic operator matrix elements: where P µ 1 = (P z , 0, 0, P z ) and P µ 2 = (P z , 0, 0, −P z ) are two momenta which approach two opposite light-like directions in the limit P z → ∞. It should be warned that this choice of normalization is not equivalent with the local matrix element P 2 ψ a Γψ b ψ c Γ ′ ψ d (0) P 1 . Γ and Γ ′ are Dirac gamma matrices, which can be chosen as Γ = Γ ′ = I, γ 5 or γ ⊥ and γ ⊥ γ 5 , so that the quark fields have leading power components on the respective light-cones. Here γ ⊥ = γ x or γ y . In principle, the combination Γ = σ µν⊥ and Γ ′ = σ µν ⊥ also gives the leading power contribution, but their matrix elements between the pion state vanish.
At large momentum transfer, the form factor factorizes through TMD factorization into TMDWFs. To motivate the factorization, one needs to consider the leading region of IR divergences in a similar way with SIDIS and Drell-Yan [50,51], and the leading reduced diagram is shown in Fig. 2. There are two collinear sub-diagrams responsible for collinear modes in + and − directions, and a soft sub-diagram responsible for soft contributions. Besides, there are two IR-free hard cores localized around (0, 0, 0, 0) and (0, b ⊥ , 0). In the covariant gauge, there are arbitrary numbers of longitudinally-polarized collinear and soft gluons that can connect to hard and collinear sub-diagrams, respectively. Based on the region decomposition, we now follow the standard procedure for the factorization [51] .
The soft divergences can be incorporated into the soft function S(b ⊥ , µ, δ + , δ − ). It resums the soft gluon radiations from fast-moving color-charges. Intuitively, soft gluons have no impact on the velocity of the fast-moving color charged partons, and the propagators of partons eikonalize to straight gauge links along their moving trajectory.
For the incoming hadron, the collinear divergences are captured by the TMDWFs for the incoming parton ψq q (x, b ⊥ , µ, δ ′ − ). However, the naive TMD-WFs contain soft divergences as well, and to avoid double counting, one must subtract the soft contribution from the bare collinear amplitude with the soft function S(b ⊥ , µ, δ + , δ ′ − ). This leads to the collinear function for the incoming direction: Similarly, for the out-going direction one obtains the collinear function Thus the explicit factorization form is conjectured as Here H F (Q 2 ,Q 2 , µ 2 ) is the hard kernel, ψ ± qq is the TMDWF, S is the TMD soft function, Q 2 = x 1 x 2 P 1 · P 2 , Q 2 =x 1x2 P 1 · P 2 . An integral over the momentum fractions x 1 ,x 2 is assumed.
Here we briefly comment on the gauge-link directions in soft functions and TMDWFs. The gauge-links along the n direction can be past-pointing. However, similar with the arguments in [52] for the SIDIS process, based on the space-time picture of collinear divergences, one can choose future-pointing gauge-links along n direction as well. With all the gauge-links being future pointing, the soft function equals to S − which is manifestly real, and the TMDWFs for the incoming and outgoing hadrons are in complex conjugation with each other. All rapidity regulators in TMDWFs and the soft functions are cancelled.
In the following, we will perform the one-loop perturbative calculation of TMDWFs, soft function, and form factor at the partonic level, and the results are presented in order.

A. TMDWFs
Since the short-distance coefficient is insensitive to the hadrons, in the calculation of TMDWFs one can replace the hadron by the partonic state. Therefore, we replace the hadron state |P by a pair of quark and anti-quark, and give the normalized definition on quark level: where the factor 1 2P + is derived from the tree-level result of the local operator matrix element, Here, the quark pair is chosen to have the same J P C with the pion, and the spin average with a Clebsch-Gordon coefficient and color average are assumed in this calculation. In Appendix A, we provide a detailed explanation of Eq. (19), and the corresponding trace formalism to derive this convention. It is necessary to mention that due to the partial conservation of axial-vector current, matrix elements in Eq. (19) are not affected by loop corrections.
According to the definition of the normalized TMDWF, one can calculate directly at the tree level, where x 0 is the momentum fraction of quark in initial state. Here, the spin average and color average are considered in this calculation. At the one-loop order, all Feynman diagrams are shown in Fig. 3. We choose the dimensional regularization d = 4 − 2ǫ to regularize the UV and IR divergences. The real diagram shown in Fig. 3(a) can be obtained as follows: The virtual diagram Fig. 3 (d) gives where L b = ln is the renormalization scale which is defined in the MS scheme. After the UV renormalization, the remanent result can be written in the form of plus function: (e) One-loop diagrams for TMDWFs. The meson state is replaced by a pair of quark and anti-quark. The first and second panels represent the real diagrams, and the third one represents the vertex diagram. The last two are virtual diagrams.
where the plus function is The summation of two Feynman diagrams exactly cancel out the infrared divergence in the delta function term.
The vertex diagram Fig. 3 (c) gives which can be rewritten as: The quark self-energy will turn the IR divergence in the second term into the UV divergence. The other two diagrams in Fig. 3 can be obtained from Eq. (23) with the exchange x ↔ 1 − x, and the details are collected in the appendix B. Summing the results in Eq. (B3, B7), and Eq. (B8), one can obtain the one loop TMDWFs where All the UV divergence term have been eliminated by composite operator renormalization. The imaginary part in Eq. (23) comes from the contribution of the Wilson line with δ − regulator, namely Fig. 3(a)(b)(d)(e), which makes opposite imaginary parts for " + " direction and " − " direction in TMDWF.

B. Soft functions and Rapidity Divergences
With a rapidity scale ζ, the rapidity divergences of the TMDWF showing in Eq. (9) can be renormalized by the on-light-cone soft functions. The soft function is defined with two on-light-cone Wilson-line cusps explicitly as where Wn is defined as Here the subscriptn give the direction of the Wilsonline, and the superscript ± in W ± n should be chosen the same as that of the WF amplitudes, and T gives the time-ordered product for quantum fields.
The soft function will be used to remove the rapidity divergence. It is interesting to notice that the soft functions can be obtained from TMDWFs with an eikonal approximation on the incoming parton lines, which re-sum the soft-gluon radiations and suffer from rapidity divergences. To ensure the scheme independence of physical TMDWFs, one needs to introduce a square root on the soft function. Since it contains two light-like directions, one can define the "physical" TMDWFs amplitudes as where y n is a dimensionless rapidity parameter for the renormalized TMDWFs. The rapidity divergences cancel between the bare TMDWFs and the soft function, which leaves a dependence of rapidity scales ζ in TMDWFs as ζ = 2(xP + ) 2 e 2yn with e 2yn = δ + /δ − . At tree-level, the matrix element in Eq. (29) only involves a N c × N c unit matrix. Therefore, it is easy to obtain S (0)± (b ⊥ , µ, δ + , δ − ) = 1.
One-loop Feynman diagrams for soft functions are shown in Fig. 4. Based on the exchange symmetry, Fig.  4 (a) and Fig. 4 (d) give the same contributions, and similar for Fig. 4 (b) and Fig. 4 (c). The results of these diagrams are given as These results in Eq. (32) and Eq. (33) do not contain the infrared divergence, because the δ + and δ − act as the infrared regulators. When q + → ∞ and q ⊥ → ∞, the soft function in Eq. (32) contains a UV divergence, which is manifested as 1/ǫ. This is similar for the kinematic region q − → ∞ and q ⊥ → ∞. The overlap of the above two kinematic regions gives the 1/ǫ 2 divergences. For the real diagrams, the results in Eq. (33) do not have UV divergence. This is due to the factor that the transverse momentum of gluon is limited by the 1/b ⊥ .
After UV renormalization, the renormalized soft functions are where S + contains an imaginary part. In the " + " direction, all Feynman diagrams in Fig. 4 contribute the imaginary part. The S − can be obtained by changing the sign in front of δ − δ + from S + . Our results are in agreement with those in the literature [53,54].
Combining the above results, we obtain the one-loop TMDWFs as whereζ = 2(xP + ) 2 e 2yn . The renormalized TMDWFs in Eq. (35) satisfies the rapidity evolution equation Substituting Eq. (35) into Eq. (36), the one-loop Collins-Soper kernel can be determined as In the above evolution equation, the CS kernel is the same with that determined from TMDDPFs.

C. Four-Quark Form Factor
In this subsection, we aim to give a complete calculation of the four-quark form factors that can be used to validate the TMD factorization scheme. At the quark level, we define the four-quark form factor as where the denominator is a normalization from two treelevel matrix elements Here, the spin average and color average are employed. Actually the spinor calculation can be implemented with a trace formalism that is described in Appendix A.
At tree level, the form factor can be directly evaluated as: It is interesting to notice that there is no UV divergence in the above equation. In addition, theses contributions are independent of the Lorentz structure Γ and the momentum fraction x 1 and x 2 . For Fig. 5 (c), the amplitude is given as: There is a three-point loop integral in this amplitude: which is rather difficult to evaluate in a brutal force way.
In the past decades, the study of mathematical properties of Feynman integrals has received increasing at-tention, and significant progress has been made in understanding the analytical behavior of multiloop Feynman integrals. One of the most powerful and advanced tools to evaluate the master integrals analytically is the method of differential equations (DEs) [55][56][57][58][59]. With the development of recent decades [60][61][62][63], this method has been widely used in various processes. Ref. [60] points out that a suitable principal integral basis (canonical basis) can be chosen in a general multiloop calculation. With the canonical basis, the corresponding DEs are greatly simplified and their iterative solutions are derived in the form of the dimensional regularization parameter. In addition, the boundary conditions of DEs can be straightforwardly determined.
To compute the integral in Eq. (44), we consider the following one-loop triangle integral family and the integral we need is G 1,1,1,0 . With the integrationby-parts (IBP) technique, one-loop QCD corrections to the real diagram of form factor are reduced into a set of integrals, named as master integrals, which are then solved using the method of DEs. The G 1,1,1,0 is evaluated as [64]: with Q ′2 = 2x 1x2 P 1 · P 2 . When the integration variable q in Eq. (45) goes to infinity, the exponential oscillation in e −iq·b ⊥ indicates the power suppression and thus there is no UV divergence in G 1,1,1,0 . The divergence in the above result is infrared.
Then the contribution from Fig. (5c) is evaluated as: Result for Fig. 5 (d) can be obtained with the replace-ment x 1 → −x 1 ,x 2 → −x 2 from Eq. (47): The vertex diagram Fig. 5 (e) gives The result for this diagram depends on the Lorentz struc-ture structure. If Γ = γ 5 or Γ = I, we obtain where Q 2 = x 1 x 2 P 1 · P 2 . Here, we use subscript '(1)' to represent the results of (pseudo) scalar structures. After absorbing the contribution from the quark self-energy one has If Γ = γ ⊥ or Γ = γ ⊥ γ 5 , the contribution is given as Here, we use subscript '(2)' to represent the results of (pseudo) vector structures. After absorbing the contribution from the quark self-energy, we have Results for Fig. 5 (f ) can be obtained from the previous results with the replacement x 1 → −x 1 and x 2 → −x 2 . If Γ = γ 5 or Γ = I we obtain Combining the above results, we obtain the complete result for the form factor with Γ = I, γ 5 If Γ = γ ⊥ , γ ⊥ γ 5 , the result is given as: A few remarks are given in order.
• The UV divergence in the I and γ 5 form factor. can be removed by the renormalization constant of scalar density operator Therefore, the renormalized form factor is • There is no UV divergence in the γ ⊥ and γ ⊥ γ 5 form factor. After some simplifications, Eq. (58) gives This is due to the fact that there is no UV divergence between the nonlocal operators, and the local ones are also free of renormalization due to the vector/axial-vector current conservation.
• An observation is that although there are infrared divergences in every diagram, summing all the results will cancel out all infrared divergences. As a result, the form factor is an infrared-safe quantity at one-loop order.

D. TMD factorization for the form factor
It has been conjectured that the form factor can be factorized into hard, collinear, and soft functions [35,37]: A rigorous proof of the factorization requests a thorough analysis of the behaviors of different dynamical modes, and in particular the cancellation of collinear and soft divergences. Though the all-order analysis is not yet given in the literature, one can use the previous results and explore the factorization at O(α s ).
At O(α s ), the form factor does not contain any infrared divergence, as shown in Eq. (58) and (60). There are infrared divergences in the TMDWFs at O(α s ), but these divergences only appear the plus function. To match the four factor, one must expand the TMDWFs and softfunctions on the right-hand side of Eq. (61) and integrate over the momentum fraction. At O(α s ), if the plus function term contributes in one of the two TMDWFs, the other quantities should take tree-level result. Thereby this term vanishes when one integrates over the momentum fraction since the hard kernel and soft function at tree-level are constants, and the other TMDWFs is a delta function. As a result, the infrared divergence on the right-hand side also vanishes. This indicates that the TMD factorization for the form-factor is valid at O(α s ).
At tree level, one has the factorization formula: from which one can obtain: In the above the arguments in H F are omitted whenever there is no confusion. It is interesting to note that tree level hard kernel can also be obtained by Fierz transformation of four-quark operators, which is detailed in Appendix A.
In a similar way, the one-loop form factor has the factorization expansion: For Γ = I or Γ = γ 5 , we have the hard kernel For Γ = γ ⊥ or Γ = γ ⊥ γ 5 , the hard kernel is calculated as The validation of TMD factorization requests a full calculation of form factors, but an interesting observation on the hard kernel can be found based on the expansion by regions. This hard function arises from the exchanges of highly offshell gluons, with a typical momentum q µ ∼ (1, 1, 1)P z . Since the amplitudes in Fig. 5 (a, b, c, d) contain an exponential factor e iq ⊥ ·b ⊥ , these amplitudes are suppressed if the momentum is hard since the factor e iq ⊥ ·b ⊥ is highly oscillating in the region Λ QCD ≪ 1/b ⊥ ≪ P z . The hard function only arises from Fig. 5 (e) and (f ), which are identical with the vertex corrections to the scalar/pseudoscalar, and vector/axialvector current. Then it can be found that the hard kernel H F (Q 2 ,Q 2 ) is determined by the spacelike Sudakov form factor as follows: where H Sud (−Q 2 ) is as given by [65].

E. Factorization Analysis based on expansion by regions
In this subsection we will adopt the expansion by regions technique and give a factorization analysis of the form factor. The analysis requests the multipole expansions of the form factor, TMDWFs, and soft function. There are a few remarks given as follows.
• Contributions from three modes are at leading power, which are hard, collinear, and soft modes with the typical momentum as: with p µ = (p + , p ⊥ , p − ), Q ∼ P z , and Λ ∼ Λ QCD .
• When the amplitudes are expanded in different regions, the lightcone divergence will show up. However, in this analysis, the δ regulator will not show up, and thus the rapidity divergence is not properly accounted for. A realistic proof in future must properly regularize the rapidity divergence.
• While the soft-function is homogeneously expanded, there are entangled contributions in TMD-WFs. Taking Fig. 3(a) as an example, this diagram contains contributions from both collinear and soft modes. When q is soft, the amplitude is simplified as: This contribution also contains a collinear contribution, and thus the Now one can perform the expansion by region for the form factor in order. For Fig. 5 (a), one can directly find this diagram is related to the vertex correction to TMDWFs, and in particular only the collinear mode contributes in this diagram: This can be derived as follows. One can make the Fierz transformation of the amplitude in Eq. (42): The structureū a (x 2 P 2 )γ ν γ 5 v d (x 2 P 2 ) is evaluated as u a (x 2 P 2 )γ − γ 5 v d (x 2 P 2 ) = 2P − 2 , and thus the γ µ in the last line of the above equation becomes γ + . Accordingly the amplitude is written as: Since at tree-level the soft function and the TMDWFs are trivial, the above amplitude takes a factorized form as Eq. (70) with ⊗ denoting the convolution over the longitudinal momentum fractions. It is similar for the conjugate diagram: in Fig. 7. One can firstly expand the amplitude with a soft gluon: When q is collinear with P 1 , one has the amplitude: . (75) This is similar with the factorization when q collinear with P 2 . Then this diagram is factorized as: where we adopted the convention for the perturbative expansion: (1/S) (1,b) = −S (1,b) . This factorization scheme is similar for the amplitude from Fig. 5 (d): . (77) For Fig. 5 (e), there are leading power contribution from the hard modes in which the full amplitude is perturbative. Contributions from the other three modes are similar with Fig. 5 (c). A sketch of the factorization is shown in Fig. 8, and the factorization formula is given as: . (78) This is similar with the factorization of Fig. 5 (f ): a) . (79) In total, the TMD factorization of the form factor at one-loop level is derived as: with the perturbative kernel up to O(α s ):

IV. QUASI-TMDWFS
The definition of quasi-TMDWFs is similar with the TMDWF, but one should replace the Lorentz structure γ + γ 5 with γ z γ 5 , and the Wilson line in quasi-TMDWF along with the z direction: where ζ z = (2xP · n z ) 2 with n µ z = (0, 0, 0, 1). Ψ ∓nz (ξ)is the field with finite Wilson line and Z E is the Wilson loop which is defined as Before presenting the result, we should mention that it is also feasible to choose the γ t γ 5 instead of γ z γ 5 in Eq. (82). But actually, there is no difference in these two Lorentz structures up to one-loop, and the verification of this behavior is given in Appendix E.

A. One-Loop Calculation
The calculation of quasi-TMDWFs requests to replace the hadron by a couple of quark and anti-quark state, At tree level, the quasi TMDWF is also a delta function, while the one-loop Feynman diagrams are shown in Fig. 9. The route C of the Wilson line W (C) = Pe ig C dsµ·A µ (s) is shown in Fig. 10. The detailed calculation of these contributions is given in Appendix D.
In the following, we take Fig. 9(b,c) as an example to illustrate the calculation. In dimensional regularization d = 4 − 2ǫ, contributions from Fig. 9 (b) and Fig. 9 (c) are given as: where L = 3 i=1 L i is the Wilson line of quasi-TMDWF: In the large L limit, we havẽ where the G is the Meijer G-function [66], and the different sign in ±i0 comes from the different Wilson line directions in Eq. (83). We should note that the third term in the above is independent of the renormalization scale. In the large P z limit, this diagram gives Summing all contributions in Eq. (D3, D4,D5) and Wilson loop in Eq. (D8), we obtain the one-loop renormalized quasi-TMDWFs: where withζ z = (2xP · n z ) 2 . The imaginary part comes from the gluon exchange between the quark field and the Wilson line, namely Fig. (9b, c).
With the above results, one can match the quasi-TMDWFs to the TMDWF as where H ± 1 ζ z ,ζ z , µ is the perturbative matching kernel. The S r is the reduced soft function defined as The S − (b ⊥ , µ, δ ± ) in the denominator is defined similar with the soft function defined in Eq. (29), but with one on-light-cone gauge-link direction alongn or n, and another off-light-cone one along n z . The reduced soft functions can also be extracted by using off-light-cone soft functions. Both the on-lightcone and off-light-cone soft functions are rapidity dependent, but the reduced soft functions are rapidity independent. The off-light-cone soft functions S ± ( b ⊥ , µ, Y, Y ′ ) are composed of two off-light-cone Wilson-line cusps. One can first define the space-like vectors asn →n Y = n − e −2Y n, n → n Y ′ = n − e −2Y ′n and the off-light-cone where the off-light-cone gauge-links Wn Y and W n ′ Y are defined as respectively. The off-light-cone soft functions are defined in a way similar to the on-light-cone soft function Eq. (29).
where Z E is Wilson loop to subtract the pinch singularities and power divergences of the off-light-cone Wilson lines. In the light-cone limit Y + Y ′ → ∞ , we have: where D is different from the on-light-cone version. Similar to the case of δ regulator, imaginary part appears in the S + case due to analyticity property. In fact, one can show that the the off-light-cone soft function depends only on the (complex) hyperbolic angle for the directions vectors from which the imaginary part can be generated. The rapidity-independent part is defined as the generalized reduced soft function: According to the properties of off-light-cone soft functions [33,67], the reduced soft function S r defined in Eq. (92) is consistent with that defined by Eq. (98). At O(α s ), the reduced soft function reads Substituting all the results known so far into Eq. (91), the matching kernel H ± 1 ζ z ,ζ z , µ is extracted as Results in Eq. (100) are in agreement with Ref. [37].

V. IMPACT ON THE REDUCED SOFT FUNCTION
In the previous section, we have validated the factorization of the form factor and determined the hard kernel in perturbation theory. A direct use of the previous results is that from Eq. (101), one can express the reduced soft function as: where the denominator term is Once the form factor and quasi-TMDWFs are simulated on the Lattice, the reduced soft functions can be determined from first-principles, and the first attempts can be found in Refs. [39][40][41][42][43][44].
In the first analyses, the tree-level result is used for the perturbative hard kernel H [40,42], while a precision de- termination requests to include the radiative corrections. Based on the lattice data on quasi-TMDWFs from Lattice Parton Collaboration with P z = 1.72 GeV, 2.15 GeV and 2.58 GeV under b ⊥ = a, 2a, 3a, 4a, 5a (a = 0.12 fm) [44], we give an estimate of the impact from O(α s ) corrections on the reduced soft function. To be more explicit, we define the ratio which directly manifests the effects on the denominator of the reduced soft function in Eq. (105). Here in H 0 , the tree-level hard kernel is used, while in H 1 , the oneloop results are incorporated. The corresponding results are shown in Fig. 11. The left panels show the results for Γ = Γ ′ = I, γ 5 , while the right ones correspond to Γ = Γ ′ = γ ⊥ , γ ⊥ γ 5 . A few remarks are given in order.
• In the lattice simulation, the Wilson line can have two directions, denoted as +L and −L. Results are consistent with each other within errors.
• Errors shown in the plots arise from the lattice data, which significantly increase with the increase of transverse separation.
• From the figure, one can see that the magnitude of QCD corrections can reach about (20 − 30)% for Γ = Γ ′ = I, γ 5 , and about −(10 − 20)% when Γ = Γ ′ = γ ⊥ , γ ⊥ γ 5 . However in the latter case, the QCD corrections to the denominator are negative, which means the corresponding reduced soft function are enhanced.
• The corrections to the denominator will decrease with the increase of P z .
• The results are insensitive to the transverse separation b ⊥ , though at large b ⊥ , the errors are too large to make a decisive conclusion.
It should be noted that the convolution in Eq. (106) involves both the longitudinal momentum fraction and the transverse separation. In general, the QCD corrections should contain the dependences on these two parameters, while results in Fig. 11 exhibit the dependences. A natural understanding of this feature is that the quasi-TMDWFs could be written as a factorized form, namelỹ Substituting this result into Eq. (106), one can see that this dependence on b ⊥ cancels in the ratio in Eq. (107). As a comparison, we adopt a phenomenological model for quasi-TMDWFΨ ′ (x, b ⊥ ) [68]: where the longitudinal and transverse distributions are entangled. We choose α = 0.197 fm, and the Gegenbauer moments a π 2 = 0.25 [69]. The behavior of R as a function of b ⊥ is shown in Fig. 12. Using the above model, we have also calculated the radiative corrections to the denominator and find that the results are about (20 − 30)% for scalar or pseudo-scalar Lorentz structure, and reach about −(10 − 20)% for vector or axial-vector structure. A dramatic difference with the results in Fig. 11 is that these results show an explicit dependence on the transverse separation.

VI. CONCLUSION
TMDPDFs and TMDWFs are important physical quantities characterizing the distributions of constituents momentum in the hadron, and reflect the nonperturbative internal structure of hadrons. In LaMET, the TMDWFs can be extracted from the first-principle simulation of a four-quark form factor and quasi distributions [35][36][37].
In the present work, a number of details are provided to understand the proposal in Refs. [35][36][37]. In particular we have explored the form factors of four kinds of four-quark operators and calculated the one-loop perturbative corrections to these quantities. In the calculation of four-quark form factors, we have adopted a modern technique based on integration by part and differential equations that can be generalized to the analysis of other nonlocal TMD quantities in LaMET. With the perturbative results, we have validated the TMD factorization of form factors and quasi-TMDWFs at one-loop level, and then extracted the O(α s ) hard function. Converting the TMDWFs to quasi-TMDWFs, the LaMET provided a "two-step" approach to the LFQ physics and achieves the goal of LFQ without performing the LFQ explicitly. Using the lattice data on quasi-TMDWFs and a phenomenological model, we have investigated the effects from the one-loop matching kernel and find that the magnitude of perturbative corrections to the soft function depend on the operator to define the form factor, but are less sensitive to the transverse separation. These results are helpful to precisely extract the soft functions and TMD wave functions from the first-principle in future.

Acknowledgment.
We thank Minhuan Chu, Yizhuang Liu, Xiangdong Ji, Kai Yan, Yibo Yang, Jialu Zhang, Jianhui Zhang and Qi-An Zhang for valuable discussions. We are grateful to Yizhuang Liu and Xiangdong Ji for carefully reading the manuscript. We thank Lattice Parton Collaboration (LPC) for allowing us to use the lattice data on quasi-TMDWFs [44]. This work is supported in part by Natural Science Foundation of China under Grants No. 12147140, No. 11735010, No. 12125503  Decay constant of a pion is parametrized by the matrix element as This provides a formal normalization for the TMDWF and quasi-TMDWFs. For the TMDWFs defined in Eq. (11), one can set b ⊥ = 0 and integrate over the momentum fraction x From this equation, one can see that the TMDWF is formally normalized to 1. Two remarks are given in order.
• Firstly, if λ = 0 and b ⊥ = 0, the Wilson line in Ψ ± n vanishes, and then the interpolating operator is local.
• Secondly, it should be noticed that the physical TMDWFs does not have to satisfy this normalization constraint. The TMDWF is valid in the hierarchy Λ QCD ≪ 1/b ⊥ ≪ 1/λ ∼ P z , however when integrating over the momentum fraction, this hierarchy is not satisfied in all kinematics region. In addition, the renormalization procedure does not commute with the integration, which also indicates that the above normalization only has a formal meaning.
At the parton level, the quark matrix element has been calculated in Eq. (39), and this can be implemented with a trace formalism. We consider a tree level matrix element: where the arrow ↑ and ↓ denote the spin +1/2 and −1/2 for quark pair. Using the spinor under the Dirac representation, one has with the coefficient c 1 = xx/2. Since this factor c 1 appears both in the evaluation of tree level and one-loop matrix elements, we can take it to be 1. Thus one can employ the matrix element 0 ψq (0) γ µ γ 5 ψ q (0) qq | tree = 1 2 sv γ µ γ 5 u ≡ 2P µ . Now we derive the tree-level matching coefficient for the form factor, based on the Fierz transformation of four-quark operators: Note that all repeated Lorentz indices in Eqs. (A6-A10) should be summed. To be complete, one should also include the color Fierz transformation: where i, j, k, l are the color indices of those quark fields ψ a , ψ b , ψ c , ψ d respectively. The second term will vanish at tree level because the index i and j are anti-symmetry for T a il . Taking Γ = Γ ′ = I as an example, we evaluate the matrix element: At tree level, there is no interaction between the quarks, and thus the four quarks can be split into two groups, each of which is related to TMDWFs: Comparing with the factorization of form factor, one can easily derive the tree-level hard kernel: This is also similar for the cases Γ = γ 5 . For Γ = Γ ′ = γ ⊥ , γ ⊥ γ 5 , notice that the form factor that we defined in the main text does not include the summation, and there is a sign difference. The tensor form factor P 2 ψ a σ µν ψ b (b) ψ c σ µν ψ d (0) P 1 seems to contribute with a leading-twist component P 2 ψ a σ µν ⊥ ψ b ψ c σ µν⊥ ψ d P 1 . However, as shown in Eq. (A10) that the Fierz transformation for tensor Lorentz structure can not generate a axial-vector Lorentz structure. Therefore, the contribution of tensor current for form factor is zero.
In addition, based on Fierz transformations one can make uses of the combinations to eliminate the power-suppressed contributions: These combinations have been used in Ref. [42].
For the vertex diagram Fig. 9 (a) we obtain: For self-energy diagram Fig. 9 (d) in large L limit we obtaiñ ψ ±(1,d) qq (x, b ⊥ , µ, ζ z ) = 1 4v (x 0 P )e −ix0P ·(λnz +b ⊥ ) γ z γ 5 dλ 2π e i(x−1)P z λ ig Then consider the Wilson loop Z E which is the denominator of the quasi-TMDWF we defined before limit. At tree-level, the operator of matrix element only involves a N c × N c unit matrix. Therefore, for the tree-level of soft function it is easy to obtain Z E (2L, b ⊥ , µ) = 1. At one-loop QCD correction as shown in Fig. 10 , the Wilson loop could be written as where the integral on the route C could be divide into four parts as C = C1 + C2 + C3 + C4 and each part itself have route C 1 (s) = −Ln z + 2Lsn z , C 2 (s) = Ln z + bs, Here the variable s ranges from 0 to 1. By computing this integral directly, we have After taking the L → ∞ limit, the final result of wilson loop can be achieved, for C i ∈ C 1,2,3,4 .
We find that the result of this diagram is independent of the Lorentz structure. The vertex diagram Fig. 9 (a) on l direction (where l = z or 0), we have the amplitude: ψ ±(1,a) qq = µ 2ǫ 0 i g 2 C F 2 (ūγ l γ 5 v) d d q (2π) d D−2 P z [(x 0 P + q) 2 q l − (x 0 P − q) 2 q l − P l q 2 ] [(x 0 P + q) 2 + iǫ][(x 0 P − q) 2 + iǫ](q 2 + iǫ) e −iq·b δ (x − x 0 )P z + q z . (E4) A brutal-force evaluation of this amplitude indicates the equivalence for the two Lorentz structures, but in the following we adopt the expansion by regions technique. Explicitly, we will demonstrate in this diagram only the collinear modes contribute.
In the quasi-TMDWFs, there are three typical models according to the decomposition of the momentum q = (q + , q ⊥ , q − ), • Hard mode with q ∼ (1, 1, 1)P z : The amplitude in Eq. (E4) contains an exponential factor e iq ⊥ ·b ⊥ , which is oscillating in the region 1/b ⊥ ≪ P z . After the integration over the q, the final result is power suppressed accordingly.
• Collinear mode with q ∼ (Q, Λ QCD , Λ 2 QCD /Q): In this region, one can find that the amplitude is O(Q), and actually the amplitudes for both structures are reduced to the TMDWF: • Soft mode q ∼ (Λ QCD , Λ QCD , Λ QCD ): In this kinematics region, one can find the power of this amplitude is O(Λ 3 QCD /Q 2 ), and namely this amplitude is suppressed.
This analysis indicates that the amplitude from Fig. 9 (a) is independent of the Lorentz structure, and moreover we have checked that the expansion by regions technique can be used to demonstrate the multiplicative factorization of the quasi-TMDWFs.