QCD Corrections of All Structure Functions in Transverse Momentum Dependent Factorization for Drell-Yan Processes

We study the one-loop correction in Transverse-Momentum-Dependent(TMD) factorization for Drell-Yan processes at small transverse momentum of the lepton pair. We adopt the so-called subtractive approach, in which one can systematically construct contributions for subtracting long-distance effects represented by diagrams. The perturbative parts are obtained after the subtraction. We find that the perturbative coefficients of all structure functions in TMD factorization at leading twist are the same. The perturbative parts can also be studied with scattering of partons instead of hadrons. In this way, the factorization of many structure functions can only be examined by studying the scattering of multi-parton states, where there are many diagrams. These diagrams have no similarities to those treated in the subtractive approach. As an example, we use existing results of one structure function responsible for Single-Spin-Asymmetry, to show that these diagrams in the scattering of multi-parton states are equivalent to those treated in the subtractive approach after using Ward identity.


Introduction
QCD factorization is an important concept for studying high energy scattering, in which both longdistance-and short-distance effects exist. With a proven factorization one can consistently separate shortdistance effects from long-distance effects. The separated short-distance effects can be safely calculated with perturbative QCD. The long-distance effects can be represented by matrix elements, which are consistently defined with QCD operators [1]. This fact allows us not only to make predictions, but also to explore the inner structure of hadrons through determining these matrix elements from experiment.
There are two types of QCD factorizations for inclusive processes. One is of collinear factorization, in which one neglects the transverse motion of partons inside hadrons at the leading twist or at the leading power. Another one is of Transverse-Momentum-Dependent(TMD) factorization, where one takes the transverse momenta of partons into account at the leading power. Using this factorization allows one to study the transverse motion of partons in hadrons and hence to obtain 3-dimensional information about the inner structure of hadrons. TMD factorization is applicable for processes in certain kinematical regions. E.g., in Drell-Yan processes, the TMD factorization can only be used if the lepton pair has a small transverse momentum q ⊥ which is much less than the invariant mass Q of the lepton pair. In this work we focus on one-loop correction in TMD factorization for Drell-Yan processes.
TMD factorization has been first studied in the case where a nearly back-to-back hadron pair is produced in e + e − -annihilations [2]. A factorization theorem in this case is established. Later, such a factorization has been established or examined in Drell-Yan processes [3], Semi-Inclusive Deeply Inelastic Scattering(SIDIS) [4,5], and has been extended to the polarized case [6]. The established TMD factorization in SIDIS and Drell-Yan processes only involves TMD quark distributions at the order of leading power q ⊥ /Q. There exist TMD gluon distributions. These distributions can be extracted from inclusive processes in hadron collision like Higgs-production [7,8], quarkonium production [9,10] and two-photon production [11]. It should be noted that studies of TMD factorization will not only help to explore the inner structure of hadrons, but also provide a framework for resummtion of large log terms with ln q ⊥ /Q in perturbative expansion with q ⊥ ≫ Λ QCD . The classical example is for Drell-Yan processes studied in [3].
Unlike parton distributions in collinear factorization at leading twist, there are many TMD parton distributions at leading twist. Structure functions, e.g., in Drell-Yan processes, are factorized with these distributions. The perturbative coefficients at tree level in TMD factorization can be easily derived. However, for reliable predictions one needs to know higher-oder corrections in the factorization. This is also important for giving physical predictions of experiments performed at different energy scales, since the dependence on the scales of perturbative coefficients appears beyond tree-level.
In Drell-Yan processes, one-loop correction of some structure functions can be obtained by studying partonic scattering and TMD parton distributions of a single parton, where one replaces each hadron with a single parton, i.e., the scattering a + b → γ * + X with a or b as a single parton state. The one-loop corrections of the studied structure functions in [6] are in fact obtained in this way. But this approach for obtaining higher-order corrections does not work for many other structure functions, e.g., the structure function for Single transverse-Spin Asymmetry(SSA) in the case that one initial hadron is transversely polarized. This structure function is factorized with the TMD parton distribution, called Sivers function [12]. If one replaces the transversely polarized hadron with a transversely polarized quark, one will always have zero results for the structure function and the Sivers function, because the chirality of a massless quark is conserved in perturbative QCD. Therefore, one needs to use multi-parton state instead of a single parton state to study those structure functions. Such a study for SSA has been done mainly in the framework of collinear factorization in [13,14,15]. The approach with multi-parton states has provided a different way to solve some discrepancies in collinear factorization of SSA [15,16].
In principle one can use these multi-parton states to study higher-order corrections in TMD factorization. Since scattering with multi-parton states is more complicated, the one-loop correction is difficult to obtain, because too many diagrams are involved. In this work we use the so-called subtractive approach to study the problem. The approach is based on diagram expansion and explained in [17,18,19]. In general, it is relative easy to find the leading order contribution to structure functions in the factorized form. At the next-to-leading or one-loop order, the diagrams, which give possible contributions, contain in general contributions from TMD parton distributions, which are of long-distance effects and need to be subtracted for obtaining the one-loop perturbative coefficients. In the subtractive approach one can systematically construct such subtractions in terms of diagrams. A comparison of the two approaches can be noticed in the following: In the approach with multi-parton states one explicitly calculates in detail the contributions of structure functions and the correspond contributions of TMD parton distributions for the subtraction. In the subtractive approach one only calculates in detail the contributions to the structure functions subtracted with the contributions of TMD parton distributions. At leading twist of TMD factorization, the symmetric part of the hadronic tensor has 24 structure functions [20,21]. With the work presented here, it turns out that the one-loop correction is the same for all structure functions. This result can be generalized beyond one-loop order.
It may be difficult to understand why the one-loop correction is the same for all structure functions. Taking SSA factorized with Sivers function as an example, the diagrams treated in the subtractive approach have no similarity to those diagrams treated with multi-parton states. We will make a comparison for a part of existing results for SSA to show that the contribution of the studied part is the same obtained from diagrams in the subtractive approach. This is in fact a consequence of Ward identity.
In Drell-Yan processes, the interpretation of the small q ⊥ is that it is partly generated with the transverse momentum of incoming partons from hadrons. In TMD factorization, as we will see in the subtractive approach, one momentum-component of partons is set to be zero as an approximation. This may give the impression that one deals here with scattering of off-shell partons and hence brings up the question if the TMD factorization is gauge-invariant. We will discuss this problem and show that the factorization is gauge-invariant.
Our work is organized as in the following: In Sect. 2. we give our notation and derive the tree-level result. In Sect. 3. we discuss the issue of gauge invariance mentioned in the above. In Sect. 4. and Sect. 5. we analyse the one-loop contributions in the factorization with the subtractive approach and give our main results. In Sect. 6. we make a comparison with a part of results derived with the subtractive approach and the existing results calculated with multi-parton states. Sect.7. is our summary. Detailed results for all factorized 24 structure functions are given in the Appendix.

Notations and Tree-Level Results
We consider the Drell-Yan process: We will use the light-cone coordinate system, in which a vector a µ is expressed as We take a light-cone coordinate system in which: h A moves in the z-direction, i.e., P + A is the large component. The lepton pair or the virtual photon carries the momentum q with q 2 = Q 2 . We will study the case with q 2 ⊥ ≪ Q 2 and that the two initial hadrons are of spin-1/2. The two hadrons are polarized. The polarization of hadron A can be described by the helicity λ A and a transverse spin vector S µ A = (0, 0, S 1 A , S 2 A ). The polarization of hadron B is described by λ B and S µ B . For convenience we also introduce two light-cone vectors: n µ = (0, 1, 0, 0) and l µ = (1, 0, 0, 0), and two transverse tensors: The relevant hadronic tensor for Drell-Yan processes is defined as: The tensor can be decomposed into various structure functions. In this work we will only give results for those structure functions which receive leading-twist contributions in TMD factorization. Taking hadrons as bound states of partons, i.e., quarks and gluons, the scattering of hadrons, hence the hadronic tensor can be represented by Feynman diagrams. Regardless how these diagrams are complicated, one can always divide each diagram into three parts: One part contains the insertion of two electromagnetic currents as indicated in Eq.(4). The other two parts are related to the hadron h A or h B . The three parts are connected with parton lines. An example is given in Fig.1a. In Fig.1a, the middle part contains the two electromagnetic vertices, the lower part is associated with h A and the upper part is associated with h B . Two quark lines connect the middle part with the part of h A and two antiquark lines connect the middle part with the part of h B . The parton lines from the part of h A only denote the contraction of Dirac-and color indices with the middle part, and the momentum flow into the middle part. The propagators associated with the parton lines are in the part of h A . The same is also for the part of h B . The middle part can be classified with the order of α s . E.g., at tree-level the middle part of   Fig.1a can be identified as: where ij stand for Dirac-and color indices. We denote the middle part in Fig.1a as H µν ij,kl (k A , k B , q), the contribution of Fig.1a as: To factorize the contribution from Fig.1a with q ⊥ ∼ Λ QCD and q ⊥ /Q ≪ 1, especially at tree level, certain approximations can be made. Because we are interested in the kinematical region of q ⊥ ∼ Λ QCD , the momenta of k µ A⊥ and k µ B⊥ can not be arbitrarily large. They are restricted as k µ A⊥ ∼ k µ B⊥ ∼ Λ QCD . Here a detailed discussion is needed to clarify what is in fact included in the hadronic parts in Fig.1b. This is also important for the comparison with the detailed calculation in Sect. 6. For this we take Γ(P A , k A ) as an example. In principle there can be the case that Γ receives contributions from large k µ A⊥ ∼ Q and also large k − A ∼ Q in Fig.1a. But these contributions can be calculated with perturbation theory. They need to be factorized out from Γ and will in general give power suppressed contributions to W µν beyond tree-level. Therefore, the dominant contributions only come from the case that Γ(P A , k A ) is only characterized with the energy scale Λ QCD , i.e., k 2 A ∼ P 2 A ∼ Λ 2 QCD . Hence one has k µ A⊥ ∼ Λ QCD and k − A ∼ Λ 2 QCD . In other word, Γ(P A , k A ) in Fig.1b is the sum of all diagrams with k µ A ∼ (1, λ 2 , λ, λ) with λ = q ⊥ /Q. Similarly, one can also find that forΓ(P B , k B ) in Fig.1b one has k µ B ∼ (λ 2 , 1, λ, λ). With the above discussion one can find the space-time picture of the hadronic matrix element in Γ. The ξ + -dependence is characterized by the small scale Λ 2 QCD . The ξ − -dependence is characterized by the scale P + A ≫ Λ QCD , and the ξ ⊥ -dependence is characterized by the scale Λ QCD . Therefore, we can first neglect the ξ + -dependence. This is equivalent to take the leading result by expanding H(k A , k B , q) in k − A . InΓ we neglect the ξ − -dependence of the hadronic matrix element. Another approximation can be made is that the leading contributions are only given by the matrix elements containing the good component of quark fields. One can always decompose a quark field as: For Γ the first term is the good component, the second term can be solved with equation of motion and gives a power-suppressed contribution. ForΓ the second term is the good component. After making these approximations, we can write the two parts as: In Eq.(8), the · · · stand for higher-twist-or power-suppressed contributions. M andM are of leadingtwist-or leading power. The quark fields in M orM are correspondingly good components. Therefore, we always have: This property will help us to extract the contributions of TMD parton distributions as we will see later.
The approximation made in the above is valid for the case that the transverse momentum q ⊥ of the lepton pair is at order Λ QCD . The correction of the approximation to the hadronic tensor is at the order q ⊥ /Q or q 2 ⊥ /Q 2 relative to the leading order. For q ⊥ ≫ Λ QCD one can make a further approximation by neglecting or expanding the ξ ⊥ -dependence in hadron matrix elements in Γ orΓ. This will lead to collinear factorizations.
We first consider the leading order given by Fig.1b. The middle part can be then given explicitly: Using the above approximated results for Γ andΓ, one obtains the hadronic tensor at leading order of α s as with q + = xP + A and q − = yP − B . According to the notation in [17], we denote the approximations made for Fig.1b to derive the above result represented by Fig.1c, it is just Fig.1b with the hooked lines. The hooked line in the lower part denotes the approximations made for Γ and the hooked line in the upper part denotes the approximations made forΓ, as indicated in Eq.(8). In the above we have worked out the contribution with a quark from h A and an antiquark from h B . Similarly, one can work out the case where a quark is from h B and an antiquark is from h A .
The density matrix M can be decomposed into various TMD parton distributions. The decomposition has been studied in [22,23,24,25]. At leading twist, the decomposition is: There are 8 TMD parton distributions at leading twist. M A is the mass of h A . In the above f ⊥ 1T and h ⊥ 1 are odd under naive time-reversal transformation. f ⊥ 1T is the Sivers function, h ⊥ 1 is called as Boer-Mulders function. For the decomposition we have implicitly assumed that the gauge links discussed in the next section is added in M.
Similarly one has the decomposition forM as: It should be noted that M andM are diagonal in colour space. With these decomposition one can work out the hadronic tensor at leading order. The results for the tensor can be represented with structure functions and each structure function is factorized with corresponding TMD parton distributions.

Gauge Invariance and Gauge Links
In the last section we have worked out the tree-level TMD factorization by considering the diagram Fig.1b. In this diagram or Fig.1a, there are only two parton lines connecting the middle part with the part of h A , and two parton lines connecting the middle part with the part of h B . From argument of power-counting, if the connection is made by more parton lines in Fig.1a or 1b, the contributions are power-suppressed but with exceptions. The exceptions are well-known. If there are gluon lines connecting the middle part with the part of h A , and those lines are for G + -gluons collinear to h A , the resulted contributions are not power-suppressed. The tree-level diagram with many collinear gluons is illustrated in Fig.2a.
It is well-known how those diagrams with many collinear gluons can be summed. The summation is achieved by introducing gauge links. In this work we follow [4] by using the gauge link along the direction u µ = (u + , u − , 0, 0) with u − ≫ u + : Finally, the summation can be represented by Fig.2b. The summation is made by re-defining: these matrices are diagonal in color-space and 4 × 4-matrices in Dirac space, and x ≥ 0. With these gauge links, TMD parton distributions will not only depend on x, k ⊥ and the renormalization scale µ but also depend on those parameters: The dependence on these parameters is controlled by Collins-Soper equation [2]. The Collins-Soper equations of the introduced TMD parton distributions can be found in [2,26]. In general one needs to add in Eq.(16) gauge links along transverse direction at infinite space-time to make density matrices gauge invariant as shown in [27]. In this work, we will take a non-singular gauge, i.e., Feynman gauge. In a non-singular gauge gauge links at infinite space-time vanish.
With the added gauge links the TMD parton distributions are gauge invariant. But there seems another problem of gauge invariance related to the tree-level result in Eq. (12). The result at the leading order seems that one can interpret the process as: One quark with the momentumk µ From the momenta one can realize that the quark and the antiquark are off-shell because of k 2 A = 0 andk 2 B = 0. One may conclude that the result is not gauge invariant because the perturbative coefficients are extracted from scattering amplitudes of off-shell partons. However it can be shown as in the following that this is not the case. We introduce two momenta which are: These are momenta of on-shell partons. Now using the properties in Eq.(10) we can derive: from this one can see that the perturbative coefficients are in fact extracted from scattering amplitudes of on-shell partons, i.e., from the annihilation amplitude q(k A )q(k B ) → γ * (q) withq µ = (q + , q − , 0, 0), indicated by the two terms in the two terms [· · ·]. The effect of transverse momenta of partons are only taken into account in the momentum conservation, i.e., in the δ-function δ 2 (k A⊥ + k B⊥ − q ⊥ ). Therefore, the tree-level result in Eq.(12) are gauge invariant. It is rather obscure to see if the tree-level result in Eq. (12) is gauge invariant from momenta carried by patrons, because the amplitudes there are constant. If we go beyond tree-level, we can see this more clearly. We consider a class of diagrams in which there is no parton crossing the cut. This case is represented by Fig.3a. After making the approximations indicated with Eq.(8) the contribution from Fig.3a can be in general written as: The contribution to perturbative coefficients from Fig.3a is obtained by subtracting the corresponding contribution of TMD parton distributions from the above contribution. Since the TMD factorization is for the leading contribution in q ⊥ /Q, one should find the leading contribution from Fig.3a. In its contribution the δ-function already gives the leading contribution. Hence we need to expand H L,R in all transverse momenta. It is easy to find the leading contribution: where · · · stand for higher order contributions in q ⊥ /Q. It is clear now in H L,R (k A ,k B ,q) the momenta of incoming partons are of on-shell. Therefore, the contribution from diagrams like that in Fig.3a is gauge invariant. This also tells us that for leading power contribution one only needs to calculate H L,R (k A ,k B ,q) with on-shell momenta. The collinear-and infrared singularities are then regularized with dimensional regularization.
We notice here that it is important and crucial to use dimensional regularization for TMD factorization here for collinear-and infrared singularities. In other word, one should set k µ A⊥ = k µ B⊥ = 0 in H L,R before performing integrations of loop momenta. One may think that one can keep a nonzero but infinite small k 2 A⊥ and k 2 B⊥ to regularize collinear-and infrared singularities. Then these singularities will appear in H L,R as terms with ln k 2 A⊥ and ln k 2 B⊥ . After subtraction of contributions from TMD parton distributions, perturbative coefficients do not contain such terms. In fact, this is not the case. This can be seen by calculating one-loop contribution of H L,R with nonzero transverse momenta of partons. The contribution will contain terms of ln k 2 A⊥ ln k 2 B⊥ . The reason for existence of such terms is that the one-loop contribution always contain double log terms. Such terms can never be subtracted with the contributions of TMD parton distributions, because the contributions of TMD parton distributions at one-loop do not have such terms. Therefore, with nonzero k 2 A⊥ and k 2 B⊥ in H L,R TMD factorization can not be made, or the TMD factorization is gauge dependent with perturbative coefficients depending on k 2 A⊥ and k 2 B⊥ . The above discussion is for diagrams without any parton crossing the cut. There are diagrams with partons crossing the cut like the one given in Fig.3b. For these diagrams there is no reason to set the transverse momenta of incoming partons to be zero as the leading approximation. These diagrams may have problems with gauge invariance. But, in TMD factorization, the contributions from diagrams like Fig.3b will be totally subtracted, as we will explicitly show at one-loop, and will not contribute to perturbative coefficients. Also, the subtracted contributions are from TMD parton distributions which are now gauge invariant, and from a gauge-invariant soft factor which will be introduced later. Therefore, there is no problem of gauge invariance with contributions from Fig.3b. Before ending this section, we briefly explain the approximation denoted by hooked lines introduced in [17], or the rule of theses lines here for TMD factorization. The approximation can be called as parton model approximation. We consider an example with a quark, which originally comes from h A and enters the annihilation into the virtual photon, as given in Fig.4. In Fig.4 the lower part is associated with h A , the upper part contains the annihilation with partons from h B . i is the Dirac index. The contribution from Fig.4 without the hooked line can be generically written as with the hooked line it means that we take the approximation for the expression as where · · · stand for contributions which will give power-suppressed contributions to the hadronic tensor and are neglected. Similarly, one can also find the rule of hooked lines for parton lines coming from the part of h B orginally.

One-Loop Real Correction
In this section we study how the real part of one-loop correction is factorized. This part corresponds to Fig.3b. We first consider Fig.5a. We will use this example to explain the subtractive approach in detail. The part needs to be subtracted from Fig.5c.
From the topology of the diagram, we can apply the parton model approximation in different places in the lower part of the diagram, e.g., one can apply the approximation for the parton lines connecting to the electromagnetic vertices directly, as given in Fig.5b, or one can also apply the approximation as given in Fig.5c and Fig.5d. It is clearly that the contribution from Fig. 5b is already included in Fig.1c, i.e., in the tree-level contribution. Therefore, in order to avoid double -counting and to obtain true one-loop contribution, one needs first to subtract the contribution from Fig.5b from Fig.5a, the remaining part will give the true one-loop correction after making the parton-model approximation for the two parton lines from the lower bubble of h A . With the approximation Fig.5a becomes Fig.5c and Fig.5b becomes Fig.5d. The true one-loop correction is then given by the contribution from Fig.5c subtracted with the contribution from Fig.5d. In the sense of the subtraction, the part of Fig.5d below the middle hooked line can be interpreted as a contribution from TMD parton distributions, because in this part all transverse momenta are at order of λ. Corresponding to the discussion for Γ after Eq.(6), this part is already included in Γ or M.
We denote the momentum carried by the gluon in Fig.5 ask A −k, and the gluon is with the polarization index ρ. k is the momentum entering the left electromagnetic vertex along the quark line from h A . The contribution from Fig.5c can be easily found as: Since we are interested in the kinematic region of q ⊥ ≪ Q, we need to find the leading contribution in the limit. For this diagram, it is easy to find the leading contribution appears if the exchange gluon is collinear to P A , i.e.: This also implies that k is collinear tok A becausek − A = 0. Using the property in Eq.(10) we obtain the leading contribution in the limit λ ≪ 1: i.e., the leading contribution has the indices µ, ν as transverse. The summed index ρ is always transverse. Again one can use the property in Eq.(10) to re-write the leading contribution in the form: where · · · stand for power-suppressed contributions which are neglected. Some integrations have been done with δ-functions. It gives q − = k − B and q + = k + . Now it is interesting to note that the part in [· · ·] is just the expression for the part in Fig.5d between the two hooked lines below the electromagnetic vertices. This part can be identified as the correction of TMD parton distributions. This fact tells that the leading contribution from Fig.5c is the same as the contribution of Fg.5d: PA PB (a) (c) Figure 6: The diagrams at one-loop for the real part.
Therefore, after the subtraction Fig.5a will not contribute to perturbative coefficients in TMD factorization at one-loop. Now we consider the contribution from Fig.6a, which also gives possible contribution at one-loop. By applying the two hooked lines in Fig.6a to parton lines nearest to the bubbles for hadrons, we obtain the contribution of Fig.6a as: herek A − k is the momentum of the exchanged gluon. It is noted that in this case the index ρ can be any of +, −, ⊥. Again we need to find the leading contributions from Fig.6a in λ → 0. These leading contributions appear in the case thatk A − k is collinear to P A or to P B , andk A − k is soft. We first consider the case that the exchanged gluon is collinear to P A , i.e., k µ ∼ O(1, λ 2 , λ, λ). The leading contribution in this case is given by: Here we use the subscript A to denote the leading contribution from the momentum region collinear to P A . It is well-known that there will be a light-cone singularity with the light-cone vector n. To avoid this we can replace the vector n with u introduced before. We can use the property in Eq.(10) to re-write the above result as: where q − = k − B and q + = k + . With the result given in the form in the above, one easily finds that this contribution is the same as that from Fig.6b. The contribution from Fig.6b is already contained in Fig.1c and should be subtracted from Fig.6a. Hence, the considered contribution in Eq.(31) is exactly subtracted and gives no contribution to the factorization at one-loop. Similarly, one also finds that the leading contribution from Fig.6a in the region where the exchanged gluon is collinear to P B is also exactly subtracted.
The interesting part is from the region where the gluon is soft, i.e., (k A − k) µ ∼ O(λ, λ, λ, λ). To analyze this part we make the substitution (k A − k) → k, i.e., now the gluon carries the momentum k. The leading contribution from this region gives: where we use the subscript S to denote the contribution from the soft gluon. We have here q + = k + A and q − = k − B . We replace the vector n with u and l with v. The color trace can be taken out by noting that the density matrices are diagonal in color space. Therefore, This contribution can be represented with Fig.6c. We note here that this contribution is not contained in Fig.1c. It can be subtracted with a soft factor as shown in [3,4]. The need of such a soft factor is also necessary for one-loop virtual correction as shown later.
If we want to subtract the soft region of the gluon from the contribution from Fig.6a, we should notice that the collinear contribution in Eq.(31) also contains the contribution from the soft region. One can easily find that the soft contribution in Eq.(31) is exactly the same as the soft contribution from Fig.6a given in Eq. (33). We can re-write the contribution from Fig.6a as: where W µν | 6aB is the contribution from the region in which the gluon is collinear to P B . This contribution contains the same soft contribution as given by Fig.6c. In the first (. . .) of Eq.(34) there are no soft-and collinear contributions. In fact the leading contribution from the first (· · ·) is zero as discussed before. In the second (· · ·) the collinear contributions are already included in Fig.1c, hence give no contribution to one-loop factorization. We define the soft factor for the subtraction of soft gluons as: At tree-level one has: At one-loop it receives contributions from Fig.7a and Fig.7b. The total one-loop correction is the sum of the two diagrams, their conjugated diagrams and those diagrams for self-energy of gauge links and for one gluon exchange between L u and L † v and between L † v and L v . The contribution from Fig.7a is: this contribution is similar to the part factorized out in Eq. (33).
With the defined soft factor we modify the factorization at tree-level in Eq. (12) as The modification will not affect the tree-level factorization. With the modification one can see that the soft contribution in the second (· · ·) of Eq. (34) is included in the soft factorS. Hence, Fig.6c will not contribute to the one-loop factorization. We can conclude that with the modification Fig.6 will not contribute to perturbative coefficients in the one-loop TMD factorization. The last diagram needs to be analyzed for the real correction is the conjugated diagram of Fig.6. The analysis is similar. The conclusion remains the same as that for the diagrams analyzed here. Therefore, we conclude that the real correction will not contribute to perturbative coefficients. This result may be extended beyond one-loop level.

One-Loop Virtual Correction and the Main Result
In this section we first analyze the correction from the virtual part and then we discuss the obtained result. For the virtual diagram there are two diagrams to be studied. One is given in Fig.8a. Another is the conjugated one of Fig.8a. After making the approximation for the lowest-and upper part, one has the contribution from Fig.8a: where we have made the parton model approximation for the two parton lines leaving the part of h A and h B . We have also taken the limit q ⊥ ≪ Q as explained for Fig.3a, i.e., the parton lines from the part of h A or h B carries the momentumk A ork B as given in Eq. (18), respectively. The gluon carries the momentum k. There are divergences in this contribution. The contribution is divergent in three regions of k: The exchanged gluon is collinear to P A or to P B , and the gluon is soft. One can work out the contributions from the three regions.
In the case that k is collinear to P A ork A , with some algebra we can write the collinear contribution in the form: where we have used the fact thatM and M are diagonal in color space and replaced the vector n with u as before. This contribution can be represented by Fig.8b. The part in [· · ·] is just the contribution from the part of the diagram in Fig.8b between the two hooked lines in the lower-half part. This contribution is already included in Fig.1c. To obtain the true correction to Fig.1c, this contribution should be subtracted. Similarly, we obtain the contribution from the region where k is collinear to P B : here the part in [· · ·] is a partM. The above contribution should be subtracted from Fig.8a too. The contribution from the region where k is small can also be worked out. It reads: This contribution is divergent and it is not subtracted by contributions of TMD parton distributions. As discussed about various contributions of Fig.6 in the last section, the collinear contributions also contain contributions from soft gluon. These soft contributions are exactly the same as given in Eq.(42).
The soft contributions will be included in the soft factor introduced in the last section in Eq.(38). The corresponding contribution to the soft factorS is from Fig.7b and reads: .
The true contribution from Fig.8a to one-loop TMD factorization can be obtained after the subtraction and the integration of k: with ρ 2 = (2u.v) 2 /(u 2 v 2 ). The · · · in the above stands for an imaginary part. The contribution with the imaginary part is cancelled when we add the contribution from the diagram which is conjugated to Fig. 8.
There is no contribution from the imaginary part in the final result of W µν at leading twist. In Eq.(44), the subtracted contributions from the [· · ·] in the first line are either already included in Fig.1c or in the soft factor. We have introduced gauge links along non-light-cone directions. This results in that in TMD parton distributions there are self-energy corrections of gauge links, corrections from gluon exchanges between L u and L † u and between L v and L † v . These corrections are all canceled by the corresponding corrections in the soft factor in Eq.(38).
From the contribution of Fig.8a, one can obtain the contribution of the conjugated diagram of Fig.8a. Summing all contributions from one-loop correction, we obtain TMD factorization of the hadronic tensor at one-loop level as: This is our main result. In the factorized form the dependence of every quantity on directions of used gauge links is explicitly indicated, and the dependence on the renormalisation scale µ is suppressed. The hadronic tensor does not depend on µ, ζ u and ζ v . From the result the one-loop correction is the same for all structure functions in TMD factorization at leading twist or leading power. The analysis here and in the last section can be generalized beyond one-loop order. The perturbative coefficient is then in fact determined by the time-like quark form factor subtracted with contributions from TMD parton distributions and from the soft factor. In our main result of the TMD factorization we have used the unsubtracted TMD parton distributions defined in Eq.(8) and the soft factorS. One can also use the subtracted TMD parton distributions and the soft factor as in [4] for the factorization without any change of the perturbative coefficient. One can also use the recently proposed definitions of TMD parton distributions in [28]. In this case, the perturbative coefficient will be changed accordingly.
With the density matrices decomposed in Eq. (13,14) one can obtain the hadronic tensor in terms of various structure functions. The results are lengthy. We give them in the Appendix, where we label the structure functions as W T U in [29]. These results agree with ours by taking ζ 2 u = ζ 2 v = ρQ 2 . The correction for the first three structure functions can be extracted by studying the partonic scattering process q +q → γ * + X by replacing each initial hadron by a single quark or antiquark. But, for remaining structure functions one can not obtain useful results by studying the partonic process. There are many reasons for it. Taking W (1) T U responsible for SSA as an example, if we replace the transversely polarized hadron h A with a transversely polarized quark and the unpolarized hadron h B with an unpolarized antiquark, one can never get nonzero result of W (1) T U of the partonic process q +q → γ * + X. The reason is the conservation of helicity in QCD.
The correction for W (1) T U in [29] is obtained through a nontrivial way. In the first step one performs a collinear factorization in the impact space with twist-3 matrix elements introduced in [30,31]. Then one adds one-loop correction. There are two relations between Sivers function f ⊥ 1T (x, k ⊥ ) and twist-3 matrix elements. The first one is to relate the second k ⊥ -moment of Sivers function with twist-3 matrix elements [32]. The second one is to relate Sivers function at large k ⊥ with twist-3 matrix elements. In this case, the relation is perturbatively determined in [33,34]. Using the second relation, one finds the one-loop correction to W (1) T U in TMD factorization. It is interesting to find that the correction to W (1) T U obtained in this way is the same as those to W (0) U U,LL,T T . From our results in Eq.(45), the correction is the same for all structure functions. At first look it may be difficult to understand this result, because the diagrams treated in [29] and the diagrams treated here are different. The difference also exists with some explicit calculations in [13,14]. We will discuss this problem in the next section.

Comparing Explicit Calculations of Multi-Parton States
In this section we take SSA as an example to discuss the problem mentioned in the above with the motivation to understand our main result. The relevant structure function is W (1) T U given in the Appendix. As explained, if we replace the transversely polarize h A with a transversely polarized quark, one will get zero results for the Sivers function f ⊥ 1T of the quark and the structure function W T U . As discussed in detail in [13], one can replace h A with a state |n = c 0 |q + c 1 |qg + · · · as a superposition of a single quark state with other states containing more than one parton. For the unpolarized h B we replace it with an antiquarkq. Then one can study the Sivers function f ⊥ 1T of the state |n and W T U of the partonic process |n +q → γ * + X. Since every quantities are of parton states, one can calculate them explicitly and directly examine factorizations.
At the leading order for the partonic process there are many diagrams [13]. One of them is shown in Fig.9a. The diagrams including Fig.9a here are already unlike the set of diagrams given in Fig.1b in the subtractive approach. In Fig.9a momenta of the partons are:p µ = (0,p − , 0, 0) corresponding to the momentum of h B , p µ = (p + , 0, 0, 0) corresponding to the momentum of h A , p µ 1 = x 0 p µ and k µ = (1 − x 0 )p µ are the momenta of the incoming quark and gluon, respectively. The short bar cutting the quark propagator means to take the absorptive part of the propagator. In fact, the short bar here is a physical cut for the absorptive part of the left amplitude.
In Fig. 9a one can nowhere to take the parton model approximation indicated by Eq.(8), or to put some hooked lines. However, in the limit of q ⊥ /Q ≪ 1, the approximation can be made. Denoting the momentum and the polarization of the gluon attached to the incoming antiquark as k g and ρ, we have for the part involving the propagator with the bar as: where γ µ is from the electromagnetic vertex. In the limit q ⊥ /Q ≪ 1, the momentum k g is at the order k µ q ∼ O(λ 2 , λ 2 , λ, λ). Hence we can approximate the above expression as: Γ ≈v(p)(ig s n ρ T a )γ µ πδ(n · k g ) ≈v(p)γ µ (ig s u ρ T a )Abs in the last step we have replaced the vector n with the vector u. Now it is clear that the contribution of Fig.9a to W T U can be factorized as given in Fig.9b with the leading order result of f ⊥ 1T of the state |n given in [13] and the leading result off 1 (y, k ⊥ ): In the limit q ⊥ /Q ≪ 1 Fig.9a is the only diagram contributing to W T U . The above discussion indicates that at the leading order the diagrams of multi-parton scattering reduce to the one treated in the subtractive approach only in the limit q ⊥ /Q ≪ 1.
At the next-to-leading order, one needs to calculate the one-loop correction of the structure function W (1) T U , and one-loop correction of Sivers function and the antiquark distribution to examine the factorization of W (1) T U . The complete calculation of these corrections will be very tedious. However, some of the one-loop real corrections of W (1) T U has been obtained in [14] for the purpose to identify the so-called soft-gluon-pole contributions in collinear factorization. One can use these results to provide certain understanding of the difference in diagrams for the derivation of our main result and for the explicit calculation presented in [14].
In Fig.9a the antiquark line does not emit any gluon which is collinear top. We consider a set of diagrams where a gluon collinear top is emitted by the antiquark lines. The contributions from this set of diagrams are of one-loop order. Some of such diagrams are given in Fig.10a,10b and 10c. There These diagrams are not similar to that given in Fig.6a. We also note here that in Fig.10a, 10b and 10c we essentially have a Glauber gluon represented by the vertical gluon line attached to a given diagram everywhere where it is possible. Because it is a Glauber gluon, Ward identity may not help.
In [14] all of the diagrams mentioned in the above has been calculated in the limit q ⊥ /Q. The leading contribution comes from the region where the gluon emitted from an antiquark line is collinear top. From [14] we have the result of this part with y = 1 after performing the integration of a loop momentum: in [14] only the divergent part with d = 4 − ǫ c is calculated. The part with y = 1 is the contribution from the region where the emitted gluon is soft. We have used the notation | y to denote the contribution from the collinear gluon. This contribution should be factorized with the leading order q T given in Eq.(48) and the one-loop real correction of the TMD parton distribution of theq given as: From the expression in TMD factorization given in Eq.(A.8) of the Appendix the discussed contribution should be factorized as: where the soft factor is at the leading order. Using the results in Eq.(48,50), one should re-produce from the right-hand side of Eq.(51) the same divergent part given in Eq.(49). It is easy to find it is indeed the case. This implies that the contributions of these more than 10 diagrams in the collinear region are factorized with the one-loop contribution off 1 (y, k ⊥ ) with y < 1. This result is in fact not strange. For this set of diagrams, Fig.10b can be easily factorized as Fig.10d. The contributions from the remaining many diagrams in the limit of q ⊥ /Q ≪ 1 can be summed via Ward identity as the contribution given by Fig.10e and Fig.10f, and hence can be factorized. Therefore, after the summation the many diagrams in the calculation with multi-parton states in the limit q ⊥ /Q ≪ 1 are equivalent to those treated in Sect. 4. This is in accordance with the discussion after Eq.(6), where it is clarified that Γ andΓ only include contributions of diagrams with transverse momenta at order of Λ QCD .
There is no calculation of one-loop virtual corrections for SSA with the multi-parton state, i.e., the virtual correction to Fig.9a. It is complicated to understand how the virtual corrections are factorized in the explicit example here. One complication comes from the one-loop correction of the left part of Fig.9a. In this part, the gluon exchanged between the gluon-and the anti-quark line is a Glauber gluon, as mentioned before. In the diagrams for one-loop correction, this Glauber gluon can be attached to several places, unlike in the case of the leading order, where the gluon can only be attached to one place as shown in Fig.9a. If one can use Ward identity for this gluon, it is rather easy to understand that the virtual corrections of the explicit example are factorized as studied in the last section. But, it is not sure if Ward identity can be used here. Besides the corrections from diagrams by adding additionally exchanged gluon in Fig.9a, there is another class of diagrams which are not generated by adding one gluon line to Fig.9a, e.g., Fig.4 in [14]. The contributions from these diagrams seem to be factorized with one-loop correction of Sivers function of the multi-parton state. To completely understand how the virtual correction in this example is factorized, a further study is needed.

Summary
By using the subtractive approach we have studied TMD factorization for Drell-Yan processes beyond tree-level. The approach is based on Feynman diagrams. In given diagrams there are nonperturbative contributions which need to be subtracted into TMD parton distributions. With the approach one can systematically construct contributions for subtracting non-perturbative effects represented by diagrams. The nonperturbative effects are only included in TMD parton distributions and the defined soft factor. After the subtraction, one obtains the perturbative contributions. We find that at one-loop the perturbative coefficients are the same for all 24 structure functions in TMD factorization at leading twist. This result may be generalized beyond the one-loop level, and the perturbative coefficient is then determined by the time-like quark form factor with the subtraction.
The QCD correction can also be obtained by studying corresponding parton processes. Replacing each initial hadron with a single parton, one can obtain the correction only for three structure functions. In this case, the diagrams are similar to those treated in the subtractive approach. For other structure functions one has to replace each initial hadron with a multi-parton state. By studying the scattering of multi-parton states, one can obtain the QCD corrections. However, the diagrams of the scattering are different than those in the subtractive approach. With existing results for the structure function responsible for SSA, one can show that the studied diagrams in the scattering of multi-parton states in the limit q ⊥ /Q ≪ 1 are equivalent to those in the subtractive approach after using Ward identity.
Our analysis can be straightforwardly generalized to the case of SIDIS. Hence, one expects that there is only one perturbative coefficient for all structure functions in SIDIS, and it is determined by the space-like quark form factor with certain subtractions. The same analysis can also be extended to TMD factorization with TMD gluon distributions in inclusive production of Higgs, quarkonium and two-photon system. One may expect that the same conclusion can be made in the case with TMD gluon distributions. Works toward this are in progress.

Appendix: The Results of the Hadronic Tensor
We organize our results of the hadronic tensor as in the following: We denote the tensor as where the index A = U, L and T denotes the contributions which do not depend on the spin of h A , the contributions depending on λ A and the contributions depending on S µ A , respectively. B denotes the similar contributions related to the spin of h B . The detailed results can be represented in terms of structure functions. We assume that the polarization of leptons in the final state is not observed. In this case, one can only measure the symmetric part of the hadronic tensor. We will only give the symmetric part of the hadronic tensor. We use the following notations: H is the perturbative coefficient given in Eq.(45). We denote the symmetric-and traceless tensor built from two transverse vectors as: In the following we give the results of structure functions in different case. For the case with AB = U U : with the following structure functions: (A.10) For AB = T L: (A.11) There are 24 structure functions in TMD factorization at leading twist. They are functions of x, y and q ⊥ . The tensor W µν U L , W µν U T can be obtained from W µν LU , W µν T U through suitable replacement. The coefficient functionsC 1,2 , C 1,2,3 and D 1,2,3,4,5 are given by: C 1 (k a , k b , q) = 1 −q 2 2q · k a q · k b − q 2 k a · k b , C 2 (k a , k b , q) = 1 −q 2 − q · k a q · k b + q 2 k a · k b , C 1 (k a , k b , q) = 1 2(q 2 ) 2 4(q · k a ) 2 q · k b − 2q 2 q · k a k a · k b − q 2 k 2 a q · k b , C 2 (k a , k b , q) = 1 2q 2 k 2 a k b · q − C 1 (k a , k b , q), D 3 (k a , k b , q) = 1 (q 2 ) 2 (q 2 ) 2 k a · k a k b · k b − q 2 (q · k a ) 2 k b · k b − q 2 (q · k b ) 2 k a · k a + (q · k a ) 2 (q · k b ) 2 , D 4 (k a , k b , q) = 1 (q 2 ) 2 − (q 2 ) 2 k a · k a k b · k b − q 2 q · k a q · k b k a · k b + 2q 2 (q · k a ) 2 k b · k b + q 2 (q · k b ) 2 k a · k a −(q · k a ) 2 (q · k b ) 2 , D 5 (k a , k b , q) = 1 (q 2 ) 2 − (q 2 ) 2 k a · k a k b · k b − q 2 q · k a q · k b k a · k b + q 2 (q · k a ) 2 k b · k b + 2q 2 (q · k b ) 2 k a · k a −(q · k a ) 2 (q · k b ) 2 , (A. 12) with k a , k b and q are transverse momenta. Here, the dot product of two vectors is defined with the metric g µν ⊥ and q 2 = q µ q ν g µν ⊥ . Through a carful comparison our results by taking the tree-level result of H = 1 agree with the existing results of factorized form factors in [20,21].