Back-to-back heavy quark pair production in Semi-inclusive DIS

The one-loop correction to heavy quark pair back-to-back production in unpolarized semi-inclusive deep inelastic scattering is given in this work in the framework of transverse momentum dependent(TMD) factorization. Both unpolarized and linearly polarized TMD gluon distribution functions are taken into account. A subtraction method based on diagram expansion is used to get finite hard coefficients. It is found the soft and collinear divergences of one-loop amplitude is proportional to tree level ones and can be expressed through several basic scalar triangle and bubble integrals. The subtraction of these divergences is spin independent. Beyond tree level an additional soft factor related to final heavy quark pair must be added into the factorization formula. This soft factor affects the azimuthal angle distribution of virtual photon in a nonperturbative way. Integrating over virtual photon azimuthal angle we construct three weighted cross sections, which depend on only three additional integrated soft factors. These weighted cross sections can be used to extract linearly polarized gluon distribution function. In addition, lepton azimuthal angle is unintegrated in this work, which provides more observables. All hard coefficients relevant to lepton and virtual photon azimuthal angle distributions are given at one-loop level.


I. INTRODUCTION
Gluon transverse momentum dependent(TMD) distribution functions take important information about the hadron structure. In high energy region, gluon as a parton usually plays more important role than quark, since the momentum fraction of parton is always very small. Thus, precise knowledge of gluon distribution functions is desired to make precise prediction of the cross section. Since gluon is spin-1, it can have different polarizations in a hadron. When the transverse momentum of the gluon in a hadron is integrated over, there is only one gluon distribution function possible for unpolarized hadron, which is the usual collinear gluon parton distribution function(PDF). In such integrated PDF the gluon is unpolarized. But when the transverse momentum of gluon is unintegrated, the transverse momentum of gluon can couple with the spin of gluon to give a description of polarized distribution of gluon in a hadron [1]. For unpolarized hadron, there exists two such transverse momentum dependent gluon distribution functions: one is denoted as G(x, p 2 ⊥ ), which reflects the distribution of unpolarized gluon in a hadron, with longitudinal momentum fraction x and transverse momentum p ⊥ , the other is H ⊥ (x, p 2 ⊥ ), which reflects the distribution of linearly polarized gluon in the hadron. So far, the second gluon distribution function, i.e., H ⊥ (x, p 2 ⊥ ) has already aroused much interest in current researches, such as the effect of H ⊥ on higgs boson production at LHC [2][3][4][5][6]. One can see e.g., [7,8] for a review on linearly polarized gluon distribution function. Various schemes are proposed to extract this function in literature. For hadron reaction, pair photon production [9], photon-quarkonium [10] or quarkonium-quarkonium [11] associated production are proposed to extract this function from a cos(2φ) azimuthal angle distribution in these reactions. For hadron reaction, the problem is the potential breaking of TMD factorization as pointed out in [12], [13]. All these proposed schemes require the detected final states are colorless(for quarkonium the heavy quark pair from hard interaction is in a color singlet). Contrast to hadron reaction, semi-inclusive deep inelastic scattering(SIDIS) for heavy quark pair or di-jet production is expect to have an exact TMD factorization, and can be used to extract some definite information about gluon TMD distribution functions. [14] have examined the effect of H ⊥ in heavy quark pair and di-jet production. For heavy quark pair production, the study in [14] is at tree level or leading order of α s with a simplified TMD formula. At this order the TMD formula just contains gluon TMD distribution functions as nonperturbative quantities. To higher order of α s , the soft radiation from final heavy quark pair will introduce an azimuthal angle dependent soft factor into TMD formula and affect the azimuthal angle distribution of virtual photon,i.e., φ q , which is used to extract H ⊥ in [14]. Because of the azimuthal angle dependent soft factor, a complete description of φ q distribution in the cross section is impossible. Instead we construct three azimuthal angle weighted cross sections, which depend on only three integrated soft factors. These weighted cross sections may help to extract H ⊥ . In [15], the TMD factorization for this process is examined at one-loop level, with H ⊥ not taken into account. It is found all collinear and soft divergences can be absorbed into gluon distribution function and a soft factor with azimuthal angle integrated. But part of the finite correction of hard coefficients is absent in [15], and the study is confined to G(x, p 2 ⊥ ). In this paper, we want to study the azimuthal angle effect introduced by the soft factor mentioned above. We will examine the factorization formula for both G(x, p ⊥ ) and H ⊥ (x, p ⊥ ) based on diagram expansion method presented in [16]. This method is different from the method in [15], which uses a single gluon to replace initial hadron. The method based on diagram expansion enables us to obtain the finite hard coefficients for various parton distributions in a systematic way. In addition, we keep lepton azimuthal angle unintegrated in this work, which can provide more observables. The structure of this paper is as follows: in Sect.II we illustrate kinematics for heavy quark pair production in SIDIS; in Sect.III, tree level factorization and resulting angular distributions are discussed, including lepton azimuthal angle distributions; in Sect.IV, the factorization formula is examined at one-loop level and the hard coefficients are calculated. In this section, the soft factors are also constructed and virtual one-loop corrections to soft factors and distribution functions are calculated; in Sect.V, the effect of soft factor is discussed and three weighted cross sections are constructed to extract H ⊥ ; Sect.VI is our summary.

II. KINEMATICS
We consider the scattering of electron and hadron where X represents undetected hadrons [17]. At order of α 2 em , this process is dominated by the exchange of a virtual photon between electron and hadron. The momentum of the virtual photon is q µ = l µ − l µ . In perturbative region, Q,Q with momenta k 1 , k 2 are heavy quark and anti-quark produced in the hard collision. In the center of mass(c.m.) frame of virtual photon γ * and initial hadron h A , we demand Q andQ are nearly back-to-back.
For convenience we define Our requirement for the final quarks then becomes where R T and K T are the transverse components of R and K, respectively. In the c.m. frame of γ * and h A , the transverse vector is relative to Z-axis. Here P A is along +Z-axis, and q is along −Z-axis. By considering only the contribution of virtual photon, the cross section we want to study can be written as dσ = 1 2s where the leptonic and hadronic tensors are In leptonic tensor we have used QED gauge invariance q µ W µν = 0 to eliminate all q µ and q ν in leptonic tensor. This will simplify our calculation. Define The phase space integration measure can be written as where ψ l is the azimuthal angle between l T and R T ; y 1,2 are the rapidities of final quark and anti-quark, respectively. Then, the differential cross section becomes where Q q is the electric charge of heavy quark in unit of the charge of electron. All information about the cross section now is contained in the hadronic tensor W µν . In next sections we will focus on the factorization of hadronic tensor.
Azimuthal angles of virtual photon and initial lepton in the transverse plane of hadron frame. γ * N frame(c.m. frame of γ * and h A ) is useful to describe the cross section, but it is not very convenient for the calculation, especially for the factorization of W µν , as we discussed later. More convenient frame is the hadron frame defined as the c.m. frame of final heavy quark and antiquark. In this frame, P A is still along +Z-axis, but virtual photon now has a small transverse momentum q ⊥ . We use light-cone coordinate representation in this paper. In this representation any vector a µ is written as (a + , a − , a µ ⊥ ), where the transverse components are relative to Z-axis and denoted by a ⊥ ; ±-components are defined by two light-like vectors n µ andn µ with n ·n = 1, i.e., a + = n · a, a − =n · a. The decomposition of relevant momenta is like The following two tensors are useful to project transverse components of a vector: and our convention for -tensor is 0123 = 1 so that 12 ⊥ = 1. For any four-vector a µ its transverse component is a µ ⊥ = g µν ⊥ a ν . Given these two frames, q ⊥ in hadron frame and K T in γ * N frame can be connected to each other. In γ * N frame, Hence, the transverse component of K T in the hadron frame is On the other hand, Combining these two equations we have This is an exact result. In the region we are considering, K T is a small quantity, then the second term with K 2 T can be ignored at leading power level.
In this work all calculations will be performed in hadron frame. To define azimuthal angles on the transverse plane we assign k 1⊥ along +X-axis, and Y -axis is defined by ⊥ -tensor, that is, Then the azimuthal angles of virtual photon and initial lepton, φ q and φ l , are defined as the angles of q ⊥ and l ⊥ relative to X-axis, respectively, as shown in Fig.1. At leading power of q ⊥ expansion, φ l is equal to ψ l , which is defined in γ * N frame. With these two vectors X µ and Y µ the transverse metric and anti-symmetric tensor can be written into another form Since n,n, X, Y form a complete basis in four-dimension space, the leptonic tensor can be expressed through these four vectors. Equivalently, one can choose P A , K, X, Y as basis. The advantage of this basis is the calculation can be done in a covariant way. The complete basis for symmetric rank-2 tensor is With this basis leptonic tensor in hadron frame is expressed as In the decomposition we have considered the constraint of QED gauge invariance, that is, P µ A and K µ must be combined intoñ µ to ensure q µ A µν i = 0. Since q · X ∼ q · Y ∼ O(q ⊥ ), QED gauge invariance is preserved at leading power of q ⊥ . One can check thatñ µ can have another representation likẽ This representation simplifies our calculation dramatically. Since hadronic tensor does not depend on φ l , this decomposition exhibits all possible lepton azimuthal angle distributions for unpolarized lepton beam.

III. TREE LEVEL AZIMUTHAL ANGLE DEPENDENCE
In hadron frame when q ⊥ is small, W µν can have a TMD factorization formula at the leading power of q ⊥ expansion. At tree level the formula reads where the gluon TMDPDF is [1] Φ αβ (x, We work in Feynman gauge ∂ µ G µ = 0 in this work. The gauge links in Φ αβ are suppressed for simplicity. There is a color summation in the definition of gluon TMDPDF. Therefore, there is a color average in the hard part. In formula eq.(21), the average factor 1/(N 2 c − 1) has been extracted, so, the hard part H µν αβ in eq.(21) contains a summation over color.
The derivation of this factorization formula for hadronic tensor has been given in [18] in detail. Under high energy limit or Q 2 → ∞, the general structure of the interaction factorizes into the form shown in Fig.2, where the two central bubbles represent hard interaction, in which all propagators are far off-shell, and the lower bubble represents the jet part of initial hadron, in which all propagators are collinear to P A . All possible soft interactions are ignored in Fig.2. These soft interactions do not appear at the leading order of α s . We will discuss them in one-loop correction. For this process, the hard scales are Q and k 1⊥ , which are taken to be the same order in this paper. Under the limit q ⊥ Q, k 1⊥ , the above factorization formula is obtained from the expansion in λ q ⊥ /Q, q ⊥ /k 1⊥ . Leading  power contribution of this process is given by collinear partons, that is, the momentum of initial gluon satisfies p µ = (p + , p − , p ⊥ ) ∼ Q(1, λ 2 , λ). According to this scaling law p ⊥ , q ⊥ are of the same order, and then the delta function δ 2 (p ⊥ + q ⊥ ) itself is already of leading power. Therefore, both p ⊥ and q ⊥ can be ignored in H µν αβ . This means H µν αβ are on-shell amplitudes. This fact ensures QED gauge invariance of hadronic tensor. For factorization in Feynman gauge gluon distributions may have a problem concerning super-leading power contribution, which appears when the gluon in Fig.2 is longitudinally polarized, i.e., G + . Such a gluon will give a 1/Λ enhancement compared to the leading power contribution we stated above. Note that the delta function in the hard part, δ 2 (p ⊥ + q ⊥ ), should not be expanded, then super-leading power contribution is given by the longitudinal gluon with p ⊥ = 0. For the one-gluon case considered here, such a contribution from longitudinally gluon vanishes due to Ward identity. But in Feynman gauge, there can be any number of longitudinal gluons connecting the central bubble and lower bubble, whose contribution is not power suppressed. For the case with two gluons connecting the central bubble and lower bubble as shown in Fig.3(a), [19] has given an explicit calculation to show the super-leading power contribution is absent even when the transverse momentum of the parton is preserved. In addition, at leading power one of the two gluons becomes gluon field strength tensor, the other is absorbed into gauge link as shown in Fig.3(b). With this conclusion we will simply take the gluon in Fig.2 as transversely polarized even at one-loop level, and will consider only the one-gluon case in our calculation. This causes no problem about the hard coefficients at least at one-loop level.
Because p ⊥ is set to zero, there is only one independent transverse momentum k 1⊥ in H µν αβ . Then previous vector basis for L µν can also be applied to H µν αβ . That is, with A µν i given by eq. (17) and Due to P-parity conservation the number of Y vector in A i B j must be even. So, there are only following 10 rather than 18 nontrivial projected hard coefficients: FIG. 4. Tree diagrams for the subprocess γ * g → QQ.
With these projected hard coefficients, the angular distributions on φ q and φ l can be obtained as where and The tree level hard coefficients can be obtained by replacing the central bubbles in Fig.2 with the diagrams in Fig.4. For the subprocess γ * + g → QQ, there are three independent variables: where p µ = xP µ A is the momentum of initial gluon. Another independent parameter we choose in our calculation is quark mass m. For convenience we defineH After some simplifications the results arẽ where From C-parity conservation H 21,22,33 are anti-symmetric, while other ones are symmetric in t 1 and u 1 . Our result satisfies this symmetry. The momentum fraction x can be obtained from in which all variables can be measured in experiment and in the last equality we have used z 1 + z 2 = 1.

IV. ONE-LOOP CORRECTION
For TMD factorization here, the relative transverse momentum of final heavy quark and antiquark is fixed. The soft divergence from virtual correction cannot be cancelled by real correction, since the phase space integration is incomplete for real correction. Usually, a soft factor is introduced to absorb such soft divergences. The operator form of the soft factor can be obtained by using eikonal approximation for soft gluons emitted by final heavy quark pair and by initial gluon, see [20] for example. The procedure is standard and the heavy quark soft factor is The definition of gauge link is with v an arbitrary vector and P the path-ordering product so that fields with smaller λ are always put on right hand side of the fields with larger λ. U v andÛ v are defined in fundamental and adjoint representations of color group, respectively. Correspondingly, G µ = G µ a T a andĜ µ = G µ aT a are gluon fields in fundamental and adjoint representations. Note thatT c ba = −if cba and T r(T a T b ) = δ ab /2 in this work. With such definitions of the gauge link our covariant derivative is D µ = ∂ µ +ig s G µ a T a . The T r[· · · ] in eq.(34) acts on matrices in fundamental representation. Then, one can check the definition eq.(34) is color gauge invariant. At order of α 0 s , S Q (b ⊥ ) is normalized to 1. The definition in eq.(34) is given in coordinate space. To get the definition in momentum space one should do a Fourier transformation, i.e., At order of α 0 s , S Q (l ⊥ ) = δ 2 (l ⊥ ). The sources of the gauge links in S Q are clear: U v1 and U v2 are obtained from the coupling of soft gluon to on-shell heavy quark or anti-quark by using eikonal approximation. Here v 1,2 = k 1,2 /m are the four-velocities of heavy quark and anti-quark, respectively;Ûṽ is extracted from the coupling of soft gluon to initial hadron or gluon. Hereṽ is collinear to the momentum of initial hadron. Ifṽ 2 = 0, S Q has a light-cone divergence, and thus is not well-defined. As a regulatorṽ is modified to be a little away from the light-cone direction but still withṽ ⊥ vanishing, i.e.,ṽ + ṽ − andṽ ⊥ = 0.
Without the soft factor the tree level TMD formula eq.(21) cannot be right. The correct one should be where H µν αβ is the hard part.S(l ⊥ ) appears in order to avoid the double counting of soft divergences, since the soft divergence in the correction to gluon TMDPDF Φ αβ is also contained in the correction to S Q . Except that now the gauge link is defined in adjoint representation,S(l ⊥ ) is the same as that defined in of SIDIS [21], i.e., In this soft factor the vectorṽ has appeared in S Q . Another vector appearing in the gauge links is As stated before, the little offshellness of v andṽ is used to regularize light-cone singularity. Other regulartors for light-cone singularity have been given by [22], [23]. The calculation procedure with these regulators is the same, so, we will not calculate the hard coefficients once more using the regulators in [22], [23].
The calculation of one-loop hard coefficients can be performed in the same way as [16]. One-loop hard coefficient H where Q ,S (1) represents the one-loop corrections to hard scattering part, gluon TMDPDF, and the two soft factors, respectively. One-loop integral in H (1) includes parton-like contribution [19], for which the loop integral is collinear to p or P A . This parton-like part has been included in tree level result and should be subtracted to avoid calculating tree level diagrams by twice. For details of subtraction one can consult [19], [16], [24]. At leading power one can show that all real corrections to hadronic tensor can be subtracted by the correction to gluon distribution and soft factors. If the final gluon is collinear to P A , [19] has shown that its coupling to heavy quarks can be absorbed into gauge link in gluon distributions. If the final gluon is soft, one can use eikonal approximation to transform the coupling of soft gluon to collinear gluon or to heavy quarks into gauge links in heavy quark soft factor S Q . The overlap between Φ αβ and S Q in soft region is subtracted by another soft factorS.
Hence, only virtual corrections contribute to one-loop hard coefficients. For virtual correction, the delta function δ 2 (p ⊥ + q ⊥ − l 1⊥ − l 2⊥ ) has been of leading power. So, p ⊥ , q ⊥ can be ignored in one-loop hard part H (1) , and then H (1) is the product of on-shell amplitudes for subprocess γ * + p → QQ. It is in this way the QED gauge invariance is preserved. If the TMD factorization formula is correct the subtracted hard coefficients H (1) f inite must be free of any infrared(IR) divergence. It should be noted that in formula eq.(37) all gluon TMDPDF, S Q andS are renormalized quantities, i.e., the UV divergences in these functions are removed by M S scheme. Thus these nonperturbative quantities all contain a renormalization scale µ. In this work, both UV and IR divergences are regularized in dimensional scheme with D = 4 − . Specially, in our scheme only loop momentum is generalized to Ddimension space, all other momenta are defined in four-dimension space. This is the four-dimensional-helicity(FDH) scheme(see [25] and references therein). This scheme is convenient for the tensor decomposition of hard part as done in eq.(23). Next we will first calculate the virtual correction to gluon TMDPDF and to the two soft factors and then present the structure of one-loop virtual correction to hadronic tensor. In the last subsection we present the explicit result of hard coefficients after subtraction.

A. Virtual correction to nonperturbative quantities
The complete gluon TMDPDF with gauge links is The virtual correction to this function is still obtained by power expansion. The diagrams contributing to virtual correction are given in Fig.5. Here we take Fig.5(a) as an example to illustrate our calculation scheme. According to collinear approximation, the leading power contribution of Fig.5(a) is from the region where the momentum of the parton going through the hooked line is collinear to P A , i.e., k µ = (k + , k − , k ⊥ ) ∼ Q(1, λ 2 , λ). Therefore, where Now in M (k) the loop integral is divergent when k ⊥ goes to zero if n = 4. But since the divergence is logarithms-like, it can by regularized in dimensional scheme. Then, M (k) is well-defined at k ⊥ = 0, and it can be expanded as Note that the hard scale in M (p + , 0, k ⊥ ) is ζ 2 = (2v · p) 2 /v 2 , so, high twist contribution is suppressed by k ⊥ /ζ. Since the delta function δ 2 (k ⊥ − p ⊥ ) is already of leading power or leading twist, only the first term in the expansion should be preserved. Thus, FIG. 6. The special vertex for the coupling of gluon field strength tensor and gauge link, Now by changing gluon field to field strength tensor the integral of above equation is just gluon TMDPDF itself. From the derivation one can see that dimensional scheme is crucial for our power expansion. Note that the nonlinear term in gluon field strength tensor also contributes. Its effect is reflected in the special vertex [26] for the coupling of gluon and gauge link, as shown in Fig.6. The total correction of Fig.5(a) is where the factor 2 represents contribution from conjugated diagrams. Note that in this expression and formulas following is IR implicitly, unless there is a special illustration. Besides, wave function renormalization for the gauge link is given by Fig.5(b) and its conjugate, the result is The sum of eq.(45) and eq.(46) is the total virtual correction to gluon TMDPDF, that is, Since the derivation does not depend on the polarization of initial hadron and the parton or gluon, the result indicates the virtual correction is the same to the two gluon TMDPDFs we consider here. Note that we will not calculate the self-energy correction to initial gluon in our following calculation for one-loop hard part, because this self-energy correction can be subtracted totally. For this reason, we do not calculate this self-energy correction to Φ αβ here.
There is no power expansion in higher order correction to soft factors, so, their virtual corrections are easy to calculate. The virtual correction to S Q (l ⊥ ) is from Fig.7. The result is v 1 and v 2 are the four-velocities of heavy quark and anti-quark, respectively. If one takes v 2 1 = v 2 2 = 1, λ becomes the usual phase space factor for final heavy quark pair, i.e., λ = 1 − 4m 2 /(k 1 + k 2 ) 2 = ρ 2 12 . The virtual correction to S(l ⊥ ) is from Fig.8, and the result is Note that in the derivation of these results we have ignored the corrections vanishing under the limit v 2 ,ṽ 2 → 0.

B. Structure of one-loop virtual correction to hard part
Next we consider the virtual correction to hadronic tensor. The central bubble in Fig.2 on the left hand side of the cut at one-loop level is given by diagrams in Fig.9, where self-energy corrections to external fermion lines are not shown, since they are trivial to be taken into account. The contribution of conjugated diagrams are not shown but taken into account in the calculation. Denoting the amplitude for the subprocess q(p) + γ * (q) → Q(k 1 ) +Q(k 2 ) as M µ α , the one-loop hard part is As an example, the tree level amplitude in Fig.4 is where α is transverse Lorentz index for initial gluon and µ is that for photon. By using this amplitude the result in eq.(31) can be obtained. From eq.(50) it is clear So, if the hard part is symmetric in (µ, α) and (ν, β), it must be real. Since the transverse momenta p ⊥ and q ⊥ are ignored in one-loop hard part, the decomposition for tree level hard part in eq.(23) can also be applied to one-loop hard part. From eq.(23), the one-loop hard part has such a symmetry automatically. So, all projected hard coefficients H ij are real. Since tree level amplitude M (0) is real, we just need the real part of one-loop amplitude. This simplifies the calculation, since absorptive part needs some caution for the i prescription in propagators in the loop. Next we discuss the IR property of one-loop amplitude. The complete amplitude is very complicated, but the IR divergent part can be shown to be simple and proportional to tree level amplitude. Here IR divergences include soft and collinear ones. These divergences can appear only in Fig.9(a,b,c,f,i). Our observation is the loop integrals with IR divergence can be expressed through three basic scalar loop integrals. Fig.9(a) contains only soft divergence, which is caused by the soft gluon exchange between the two quarks. Such soft divergence is contained in following integral, FIG. 9. One-loop virtual corrections to the amplitude in hard scattering part, where the self-energy corrections to external heavy quark and anti-quark is not shown, since they are trivial to be taken into account.
The soft divergence of such box integral can be expressed through scalar triangle integrals [27]. In soft region with l µ = (l + , l − , l µ ⊥ ) ∼ Q(λ, λ, λ), (l + p − k 2 ) 2 − m 2 is offshell. Thus setting l = 0 in this propagator does not affect the IR divergence. This off-shell propagator can be decomposed as The second term gives IR finite contribution. Then, All soft divergence of the box integral now is contained in the triangle integral CTri1(s 1 ). Note that CTri1(s 1 ) is symmetric in k 1 , k 2 or t 1 , u 1 . By exchanging k 1 and k 2 , the soft divergence of Fig.9(b) is obtained. After some simplifications the sum of the divergent parts of Fig.9(a,b) is The quantity in {· · · } is rightly the tree level amplitude. Hence, For Fig.9(c), both soft and collinear divergences are contained in the loop integral, and the overlap of these two divergences makes the extraction of IR part nontrivial. The divergent loop integrals for this diagram are with In I Boxc , collinear divergence comes from the region where l µ is collinear to p µ , while soft divergence comes from two regions: one is l µ is soft, the other is (p − l) µ is soft. In I (1,2) Boxc , collinear divergences come from the same region as Boxc , but the soft divergence only comes from the region where (p − l) µ is soft. Note that the integral with (l + ) 3 in the integrand is absent in the amplitude, although this integral is also IR divergent.
As we will show, the divergent part of Box-c can still be expressed through triangle integrals. First, setting l = 0 in D 2 gives where . = means the equality holds up to IR finite corrections. Now l · k 2 = l + k − 2 + l − k + 2 + l ⊥ · k 2 , but l − and l ⊥ are suppressed by λ or λ 2 in soft or collinear regions, so, l − and l ⊥ can be dropped. Then where we have used q µ = −p µ + k µ 1 + k µ 2 . Now D 0,1 in the numerator above can be dropped, because they are suppressed by λ or λ 2 in soft or collinear regions. Then, substituting this result back to eq.(60) one gets In this way we express the IR part of Box-c through two triangle integrals. As expected, this expression is symmetric in k 1 , k 2 . In the above we have also defined the finite part of this box integral as DBox2[t 1 , u 1 ]. Similarly, one can show With the same method, the extraction of IR divergence from the triangle diagrams Fig.9(f) becomes easy, and the results are The integrals with (l + ) 2 or (l + ) 3 do not appear in the amplitude, although they are divergent. Other integrals are IR finite for this triangle diagram. By exchanging k 1 , k 2 the divergent part of Fig.9(i) can be obtained. With the obtained IR divergence, the amplitude from Fig.9(a,b,c,f,i) can be expressed through four basic scalar integrals J 01 , J 012 , J 013 and CTri1(s 1 ). After a lengthy calculation the result is put into a very neat form, which is proportional to tree level amplitude! Besides, the self-energy corrections to external fermion lines also contain IR divergence. Then, the complete IR divergent part of the amplitude is with This is one of our main results. The derivation of this result is independent of the polarizations of external gluon, thus this result indicates the IR divergence, including soft and collinear ones, is spin independent. According to eq.(50), the IR divergence of hard part is where the factor 2 indicates the contribution from conjugated diagrams.
Besides IR divergence, the hard part also contains UV divergence, which is subtracted by the counter terms of lagrangian. After the subtraction of UV divergence, the hard part gives a µ dependence, with µ the renormalization scale. We will separate the µ dependence in the following.
Notice that the electro-magnetic current j µ in hadronic tensor is conserved, so, the UV divergence from vertex correction is cancelled by self-energy correction to fermions. This means the UV divergences of Fig.9(e,h) are cancelled by Fig.9(j,k). So, the sum of these four diagrams does not generate µ dependence. Apart from the µ dependence in wave function renormalization constant for fermions, i.e., Z 2 , the remaining µ dependence can only come from Fig.9(d,f,g,i). The UV divergence of these diagrams reads With the µ dependence from Z 2 , the total µ-dependence of the amplitude is where · · · represents µ independent part. With renormalization scale separated explicitly, the hard part can be written as with where we have chosen the reference scale of µ 2 as Q 2 so that H 3 is definite. These three parts make the structure of one-loop virtual correction to hadronic tensor clear. At last, it should be pointed out that only in pole mass scheme, the µ dependence of Fig.9(j,k) is proportional to tree level amplitude. In pole mass scheme the mass counter term in lagrangian is chosen so that renormalized fermion self-energy Σ R (/ p, m) = 0 when / p = m. Then the renormalized self-energy is written as with From the structure of self-energy correction, ln µ part is proportional to tree diagrams. As stated before this µ dependence is cancelled by that of Fig.9(e,h). This is an advantage of pole mass scheme.

C. Subtraction and finite hard coefficients
According to our subtraction scheme in eq.(39) and results in eq. (47,48,49,67), the finite hard coefficient is given by where we have used the fact and the ln µ 2 is from Z 2 and J 01 . The explicit expression of IR divergent scalar loop integrals are calculated by standard Feynman parametrization. One can also find general expressions in e.g. [27], [28] and then continue the result into DIS region. Generally, the result in DIS region is very simple. For completeness we list the result of these integrals here as and where and In DIS region, y 13 , y 23 < 0 and y 12 , r b > 0. In these integrals we have ignored all absorptive parts. With these integrals, the explicit W IR can be obtained as The expression of W Φ has been given by eq.(47), then we have It can be seen the IR divergence from hard part cannot be subtracted by the correction to gluon TMDPDF alone. The TMD formula without soft factors must be wrong. At one-loop level, the two soft factors have been calculated in eq.(48) and eq.(49). We have Now it is very clear that the difference between W IR − W Φ and W SQ + W S is IR finite, i.e., In this result, the UV counter terms for gluon TMDPDF and soft factors have been added and the scale µ is for UV renormalization. In eq.(83) the first line in the equality is given by wave function renormalization, that is, After UV divergences are removed the wave function renormalization constants are So, Now, substituting eq.(83) into eq.(74), the finite hard coefficient is Moreover, our decomposition for hard coefficient in eq.(23) is effective to all orders of α s if only virtual correction is involved. So, H µν 2αβ can still be projected into 10 scalar hard coefficients H ij 2 as we done for H µν(0) αβ . There will be no other new tensor structure in higher order correction. That is, with the projection in eq.(23). There is not any new tensor structure appearing in higher order correction, since only virtual correction contributes. Eq.(87,88) is one of our main results. From the derivation, W IR and W Φ are independent of the polarizations of parton or gluon, so, W f is the same to all gluon TMDPDFs, including those defined by polarized hadron as given in [1]. H 2 depends on the polarizations of gluon, but it is IR finite and µ independent. The finiteness of H 2 is justified explicitly according to eq.(71). In our result, it is expressed through some basic IR finite scalar loop integrals, including box integrals DBox1, 2, triangle integrals CTri3, 4, 5 and various bubble integrals b 0 , B v and B m . These integrals are illustrated in Appendix.. In our result, the projected hard coefficients H ij 2 according to eq.(23) are given in a mathematica file subted2.m. These results are very lengthy and cannot be shown here. In subted2.m these projected hard coefficients are stored as a list with normalization factors All of these results are expressed through the invariants s 1 , t 1 , u 1 and m appearing in tree level result eq.(31). There are two color factors in the stored result, i.e., In the final part it is interesting to point out that our factorization formula eq.(37) associated with the finite hard part given by eq.(87) is µ independent to order of α 2 s . This can be illustrated in the following. Because real correction is UV finite, µ dependence is purely generated by virtual corrections. Then, the µ dependence of S Q (l ⊥ ) andS(l ⊥ ) can be obtained from eq.(48) and eq.(49), respectively. They are For gluon TMDPDF, additional self-energy correction to gluon with momentum p should be added into virtual correction. Of course, this correction does not affect the finite hard part, because it cancels the same gluon self-energy insertion in hard scattering for hadronic tensor. The wave function renormalization constant in Feynman gauge is in M S scheme. Then, from eq.(47) we have Recall that tree level hard coefficients are proportional to α s , which has a µ dependence FIG. 10. Azimuthal angles φ b and φq in hadron frame.
Now the hadronic tensor to O(α 2 s ) can be written as where · · · part is O(α s ) and µ independent; ⊗ means the transverse momentum convolution. From this equation and evolution equations given above, one can check quickly So, the factorization formula is µ independent to order α 2 s . Notice that here we have used the fact that pole mass of heavy quark is µ independent.

V. AZIMUTHAL ANGLE DEPENDENCE
Now we have checked the corrected TMD formula for hadronic tensor, i.e., eq.(37) at one-loop level, and the finite hard coefficients are given explicitly. One may think that the azimuthal angle dependence about φ q can be obtained by following the same procedure as tree level formula. But unfortunately, S Q depends on v 1⊥ , which will change the tree level φ q dependence in a nonperturbative way after integration over the transverse momenta in eq.(37). Physically, this change of φ q dependence is caused by the multiple emissions of soft gluons from final heavy quark pair. To proceed we have to integrate over φ q , and it is best to do this in coordinate space. DefineΦ Note that b 2 ⊥ = b ⊥ · b ⊥ < 0. The formula eq.(37) becomes Remember that in hadron frame, i.e., the rest frame of final heavy quark pair, v 1⊥ is along +X-axis, then all azimuthal angles are relative to v 1⊥ , as shown in Fig.10.
Since L µν does not depend on φ q , it is convenient to integrate over φ q in W µν rather than the cross section. Notice that due to scaling invariance of That is,S Q is a function of cos 2φ b . So,S Q is unchanged under the transformation φ b → φ b + π. Due to this fact we have This property ofS Q is valuable in our following analysis. We choose to study three quantities with φ q integrated: are the basic tensors we introduced in the decomposition of leptonic tensor, i.e., eq. (17). For simplicity, the quantity concerning H ⊥ is reorganized intoH Due to the integration over φ q , φ b now can be integrated. As a result, three integrated soft factors are involved, i.e., Generally, the quantities weighted by cos(nφ q ), sin(nφ q ) can be constructed, but more φ b integrated soft factors will be introduced. This may not help to extract H ⊥ . From these weighted hadronic tensors, the corresponding cross sections can be obtained. Define weighted cross section as σ cos nφ q = 2π 0 dφ q cos(nφ q ) dσ dψ l dydx B dy 1 dy 2 d 2 K T d 2 R T .
We have Here H ij = H (0)ij + H (1)ij f inite , which are given in eq.(31,88). These three weighted cross sections are our main results. We can see clearly that σ 1 also depends on H ⊥ even if lepton azimuthal angle φ l is integrated over, due to nonzeroS (2) Q . This feature is unexpected from tree level TMD formula eq. (21). In addition, the soft factorS Q is process independent and can be absorbed into gluon distributionG orH ⊥ , which results in so called subtracted TMDPDF [21]. ButS Q is process dependent and cannot be eliminated. Obviously,S Q is similar to fragmentation functions in SIDIS. The difference is fragmentation functions can be extracted from e + e − experiment, butS Q cannot be extracted in this way, because it contains a gauge link related to initial gluon. The detailed knowledge ofS Q is beyond the scope of this paper, and we will try to study its effect on resummation in future.

VI. SUMMARY
In this paper we first derive the angular distributions of heavy quark pair back-to-back production in SIDIS based on tree level TMD formula. Then we examine the one-loop correction to this formula. At one-loop level a special soft factor for final heavy quarks should be complemented. From our previous studies we have known that real correction does not contribute to higher order hard coefficients. Therefore, we calculate only virtual corrections to the hadronic tensor and various nonperturbative quantities in factorization formula. Really we find the IR divergence in hadronic tensor can be absorbed by these nonperturbative quantities. As a result, we give explicit form of finite hard part, including renormalization scale dependence. Interestingly we find the IR divergent part of one-loop amplitude for heavy quark pair production is proportional to the tree level amplitude and can be expressed through standard triangle and bubble loop integrals. This feature ensures the subtraction in polarized scatterings can also be done in the same formalism. In order to present the azimuthal angle dependence about virtual photon, i.e., φ q in hadron frame, we project the hard part into ten scalar hard coefficients. However, at one-loop level, the appearance of soft factor S Q affects the azimuthal angle dependence in a nonperturbative way, which makes the explicit φ q dependence at cross section level become unclear. In order to extract some information about gluon TMDPDFs, especially about linearly polarized gluon distribution H ⊥ , we construct three φ q -weighted cross sections σ 1 , σ cos 2φ q , σ sin 2φ q , which depend on only three integrated heavy quark soft factorsS (0),(2),(4) Q . We expect future experiments such as EIC [29] can give some constraints on these three soft factors and gluon distribution H ⊥ . For TMD factorization, an important issue is the resummation of relative transverse momentum of heavy quark pair in this process. But due to the many hard scales in this process, e.g., Q 2 , R 2 ⊥ , and the scales introduced by soft factors like ζ 2 , 4v 1 ·ṽv 2 ·ṽ/ṽ 2 , the resummation will be nontrivial. We want to study this issue in another paper. In appendix we present all involved IR finite loop integrals, which are used to express the finite hard part H ij 2 .