Factorization and its Breaking in Dijet Single Transverse Spin Asymmetries in $pp$ Collisions

We study factorization in single transverse spin asymmetries for dijet production in proton-proton collisions, by considering soft gluon radiation at one-loop order. We show that the associated transverse momentum dependent (TMD) factorization is valid at the leading logarithmic level. At next-to-leading-logarithmic (NLL) accuracy, however, we find that soft gluon radiation generates terms in the single transverse spin dependent cross section that differ from those known for the unpolarized case. As a consequence, these terms cannot be organized in terms of a spin independent soft factor in the factorization formula. We present leading logarithmic predictions for the single transverse spin dijet asymmetry for $pp$ collisions at RHIC, based on quark Sivers functions constrained by semi-inclusive deep inelastic scattering data. We hope that our results will contribute to a better understanding of TMD factorization breaking effects at NLL accuracy and beyond.


I. INTRODUCTION
There has been a strong interest in correlated dijet production in various hadronic collisions [1][2][3][4][5][6][7][8][9], where the two jets are produced mainly in the back-to-back configuration in the transverse plane, A + B → Jet 1 + Jet 2 + X . (1) Here A and B represent the two incoming hadrons with momenta P A and P B , respectively. The azimuthal angle between the two jets is defined as φ = φ 1 − φ 2 with φ 1,2 being the azimuthal angles of the two jets. In the leading order naive parton picture, the Born diagram yields a delta function at φ = π. One-loop gluon radiation will lead to a singular distribution around φ = π. This divergence arises when the total transverse momentum of the dijet (imbalance) is much smaller than the individual jet momentum, q ⊥ = | k 1⊥ + k 2⊥ | ≪ |k 1⊥ | ∼ |k 2⊥ | ∼ P T , where large logarithms appear at every order of the perturbative calculation.
In the kinematic region q ⊥ ≪ P T , the appropriate resummation method that needs to be applied is the so-called transverse momentum dependent (TMD) resummation or the Collins-Soper-Sterman (CSS) resummation [10]. There have been several theoretical efforts to resum the large logarithms for this process [11][12][13][14][15][16][17][18]. The differential cross section can be written as, where dΩ = dy 1 dy 2 dP 2 T represents the phase space of dijet production. Here y 1 and y 2 are the rapidities of the two jets, P T is the leading jet transverse momentum, and q ⊥ the imbalance transverse momentum between the two jets as defined above. Moreover, σ 0 is the overall normalization of the differential cross section. The first term on the right hand side, W ab→cd , contains the all order resummation and the second term, Y ab→cd , takes into account fixed order corrections. At next-to-leading logarithmic (NLL) order, the resummation for W was conjectured to take the following form [16,17] W ab→cd (x 1 , x 2 for each partonic channel ab → cd, where Q 2 =ŝ = x 1 x 2 S, representing the hard momentum scale. In addition, we have b 0 = 2e −γ E , with γ E being the Euler constant. The f a,b (x, µ) are the parton distributions for the incoming partons a, b, and x 1,2 = P T (e ±y 1 + e ±y 2 ) / √ S are the fractions of the incoming hadrons' momenta carried by the partons. In the above equation, the hard and soft factors H and S are expressed as matrices in the color space of the partonic channel ab → cd, and γ s ab→cd are the associated anomalous dimensions for the soft factor. The Sudakov form factor S Sud resums the leading double logarithms and the universal sub-leading logarithms, where R 1,2 are the jet radii of the two jets, respectively. In practice the jets are of course reconstructed with the same radius R but to clarify the structure of our calculation we use two different radii R 1,2 to differentiate between the dijets. At one-loop order, A = C A αs π , B = −2C A β 0 αs π for a gluon-gluon initial state, A = C F αs π , B = −3C F 2 αs π for a quark-quark initial state, and A = (C F +C A ) 2 αs π , B = ( −3C F 4 − C A β 0 ) αs π for a gluon-quark initial state. In addition, D 1,2 = C A αs 2π for a gluon jet and D 1,2 = C F αs 2π for a quark jet, respectively. Here, β 0 = (11 − 2N f /3)/12, with N f being the number of effective light quarks.
The resummation formula in Eq. (3) was obtained in Refs. [16,17] by a detailed analysis of the soft gluon radiation at one-loop order. The leading contributions from soft gluon radiation can be factorized into the associated TMD parton distributions and can be resummed by solving the relevant evolution equations. At NLL, the soft gluon radiation is factorized into the soft factor S which is given by a matrix in the color space of the partonic channels. The matrix form of the factorization is the same as was found for threshold resummation for the dijet production in proton-proton collisions [19][20][21][22][23][24].
It is known that TMD factorization in dijet production in hadronic collisions is highly nontrivial and that there are potential factorization breaking effects [25][26][27][28][29][30][31][32][33][34][35]. First, nonglobal logarithms (NGLs) [36,37] start to contribute to the cross section at two-loop order. It has been shown that they cannot easily be included into a factorization formula, although numerical simulations can be made and their contribution can be taken into account [11,12]. In addition, TMD factorization will be explicitly broken at three-loop order for the unpolarized cross section. This leads to a modification of the coefficient A (3) in the above Sudakov form factor [17, 27, 32-35, 38, 39].
Factorization breaking effects are particularly evident for the single transverse spin asymmetry (SSA) in dijet production [27], ∆σ(S ⊥ ) = (σ(S ⊥ )−σ(S ⊥ ))/2, where S ⊥ represents the transverse polarization vector for one of the incoming nucleons. The SSA for this process is expressed as ∆σ(S ⊥ ) ∝ ǫ αβ S α ⊥ q β ⊥ , i.e., the total transverse momentum of the two jets q ⊥ will have a preferred direction [25]. This asymmetry is sensitive to the parton's Sivers function where the transverse momentum distribution is correlated with the transverse polarization vector [40]. In Refs. [29,31], all initial/final state interaction contributions to the SSA were factorized into a complicated gauge link structure associated with the quark Sivers function for the polarized nucleon. However, for the double spin asymmetries involving two Sivers functions, it was shown explicitly that the generalized gauge-link approach to TMD factorization does not apply [28].
Ref. [26] provided an understanding of the SSA from the twist-three framework where the Qiu-Sterman-Efremov-Tereyav matrix elements are the basic ingredients [41][42][43][44][45], and the high momentum Sivers function is generated by collinear gluon radiation. In particular, it was shown in Refs. [26,42] that the collinear gluon radiations parallel to the incoming hadrons can be factorized into the associated TMD parton distribution functions. It was also suggested that a factorization formula similar to that in the unpolarized case may hold for the single spin dependent differential cross section [26], × S ab→cd (λ ⊥ ) H Sivers ab→cd (P 2 ⊥ ) c δ (2) ( p 1⊥ + p 2⊥ + λ ⊥ − q ⊥ ) .
Here f

⊥(SIDIS)
1T b and f (SIDIS) a denote the transverse-spin dependent TMD quark Sivers function and the unpolarized TMD parton distribution, respectively. These TMD parton distribution functions were chosen following their definitions in semi-inclusive deep-inelastic scattering (SIDIS) with future pointing gauge links. Although it was not explicitly shown in Ref. [26], a matrix form of the factorization was suggested, where H Sivers ab→cd and S ab→cd are partonic hard and soft factors and the [ ] c represents a trace in color space between the hard and soft factors, similar to the unpolarized case in Eq. (3).
In order to check the factorization formula of Eq. (5), it is important to carry out the calculation of soft gluon radiation. Soft gluon emissions contribute in a nontrivial way to the factorization formula. In particular, it will be crucial to show that these contributions can be included in the soft factor in the matrix form of the factorization formula. The goal of the current paper is to derive the soft gluon radiation contribution at one-loop order. As mentioned above, in Ref. [26] it was shown that collinear gluon radiation associated with the incoming nucleons can be treated following the general factorization arguments. This indicates that factorization holds in the leading logarithmic approximation (LLA). However, in order to obtain also all the subleading logarithmic contributions, we need to consider the soft gluon radiation as well. After including soft gluon radiation, we will obtain the complete double logarithmic result.
Our calculations presented in this work show that the factorization and resummation is expected to be valid at LLA. However, factorization breaking effects will emerge at NLL accuracy, in the sense that the contributions from soft gluon radiation cannot be factorized into the same soft factor as for the unpolarized case. This implies that beyond the LLA, we do not have a factorization formula for ∆σ as in Eq. (3), at least not in the standard way with a spin independent soft factor.
In the LLA, we can express the spin dependent differential cross section in terms of the Fourier transform b ⊥ variable, Here, we neglect the Y -term contribution compared to the unpolarized case above. In this work we show that the leading logarithmic factorization of W T β takes the form where S T Sud (Q 2 , b ⊥ ) can be written in analogy to Eq. (4), Here A, B, D 1,2 are the same as in Eq. (4). In the above equation, T F is the Qiu-Sterman matrix element which is also related to the transverse momentum-moment of the quark Sivers function. It is defined as follows where F µν represents the gluon field strength tensor. From the leading order derivation, we have 1 . For our soft gluon calculation at leading power, we take the correlation limit, i.e., we neglect all power corrections of the form q ⊥ /P T . In this limit, the leading double logarithm is proportional 1/q 2 ⊥ × ln(P 2 T /q 2 ⊥ ). We will show that these contributions will be consistent with the resummation formula of Eq. (7). Some sub-leading logarithmic terms can be factorized in this form as well. These include collinear gluon radiation associated with the incoming hadrons and the final state jets. The former can be resummed by including the scale evolution of the integrated parton distribution and the Qiu-Sterman matrix elements, e.g., by evaluating these distributions at the scale The latter are taken into account by the ln(1/R 2 ) terms in the Sudakov form factor S Sud of Eq. (4). Therefore, Eq. (7) is an improvement of the leading logarithmic approximation, to which we will refer as LLA ′ in the following.
The remainder of this paper is organized as follows. In Sec. II, we will briefly review the soft gluon contribution to unpolarized dijet production. The basic formalism, including the Eikonal approximation, the phase space integrals to obtain the leading contributions, and the subtraction method to derive the soft gluon radiation associated with the final state jets will be introduced. Sec. III contains the main new derivations of this work. We will carry out the calculation of the soft gluon radiation for the spin dependent differential cross sections. We introduce the general framework, the twist-three collinear expansion, and derive the soft gluon radiation amplitude in this formalism. We apply these techniques to different partonic channels and demonstrate that the leading double logarithmic contributions factorize, and we verify our resummation formula at LLA ′ . In Sec. IV, we consider an example and show that factorization breaking effects appear at NLL, in particular, from soft gluon radiation that does not belong to the incoming hadrons and final state jets. In Sec. V, we will present phenomenological results for the single spin asymmetries in dijet production at RHIC, and compare to recent STAR data [2,46]. Finally, we will summarize our paper in Sec. VI.

II. BRIEF REVIEW OF SOFT GLUON RADIATION FOR THE UNPOLARIZED CASE
Dijet production at the leading order can be calculated from partonic 2 → 2 processes, where p 1,2 and k 1,2 are the momenta of the incoming and outgoing two partons, respectively. Their contributions to the cross section can be written as where the overall normalization of the differential cross section is σ 0 = α 2 s π s 2 . The partonic cross sections h (0) for all the production channels depend on the kinematic variablesŝ = (p 1 + p 2 ) 2 ,t = (p 1 − k 1 ) 2 andû = (p 1 − k 2 ) 2 . As mentioned above, at the leading order, they contribute to a delta function setting q ⊥ = 0, which corresponds to the back-to-back configuration of the two jets in the transverse plane. For soft gluon emissions, we can apply the leading power expansion and derive the dominant contribution from the Eikonal . 1: Soft gluon radiation contributions to the finite imbalance transverse momentum q ⊥ : (a) initial state radiation, and (b), (c) final state radiation. Since we have chosen the gluon polarization vector along p 2 , there is no gluon radiation from the line with momentum p 2 .
approximation [15,17]. For example, for the outgoing quark, antiquark and gluon lines, we obtain the following factors in the Eikonal approximation: respectively, at one-loop order. Here g is the strong coupling and the k i represent the momenta of the outgoing particles. For incoming quark, antiquark and gluon lines, we have, respectively, where p 1 corresponds to the momentum of the incoming particle. Following Ref. [15,17], we choose physical polarizations of the soft gluon along the incoming particle with momentum p 2 , so that the polarization tensor of the radiated gluon takes the following form: This choice will simplify the derivation since there is no soft gluon radiation from the incoming parton line p 2 . Therefore, as shown in Fig. 1, the leading contributions come from the initial state radiation from the line with momentum p 1 , and the final state emissions from the lines k 1 and k 2 . The contributions by these diagrams can be evaluated by taking the amplitudes squared of the Eikonal vertex with the polarization vector of the radiated gluon contracted with the above tensor. This leads to the following expressions for the soft gluon radiation contributions, Here S g (p, q) is a short-hand notation for When we integrate out the phase space of the radiated gluon to obtain the finite transverse momentum for the dijet imbalance, we have to exclude the contributions that belong to the jets. Therefore, only gluon radiation outside of the jets with radius R 1,2 contributes. These diagrams have been calculated in Refs. [17], where an offshellness was considered to regulate the collinear divergence associated with the jet within the narrow jet approximation [47,48]. In Ref. [49,50], a subtraction method was employed to derive the soft gluon radiation contribution. For completeness, we show details of the derivation in the Appendix. Here, we list the final results: Compared to the results in Ref. [17], the above results differ by a term proportional to ǫ π 2 /6. This is a result of the approximation made in Ref. [17].

FIG. 2:
Leading order diagrams for the initial and final state contributions to the SSA in dijet production. The red bars in these diagrams indicate the propagators that produce the necessary phase for the SSA.

III. SOFT GLUON RADIATION FOR SSA: LEADING LOGARITHMIC CON-TRIBUTIONS
In this section, we will investigate the soft gluon radiation contribution to the SSA in dijet production. The leading order analysis and collinear gluon radiation contributions have been studied in Ref. [26]. In the following, we will first review the leading order results and then derive the soft gluon radiation contribution.

A. Leading Order Results
The leading order results of Ref. [26] can be transformed into the factorization formula of Eq. (7). For convenience, we show the diagrams that contribute to the SSA from initial and final state interaction effects in Figs. 2. As demonstrated here, the SSA phases only come from the gluon attachments to the initial/final state partons. The leading order results derived in Ref. [26] can be obtained from Eq. (5) by setting S ab→cd (λ ⊥ )H Sivers ab→cd (P 2 T ) c ≡ H Sivers ab→cd (P 2 T ). After taking the Fourier transform to impact parameter b ⊥ -space, we find the following leading order result: The hard part is written as where i labels the different contributions to the hard factors h i by the various Feynman diagrams. Here the factors C i I are for the initial state interaction for the single-spin dependent cross section, and C i F 1 and C i F 2 are for the final state interactions when gluon is attached to the lines with momentum k 1 and k 2 , respectively. The explicit expressions for and h i are given in Ref. [26].
Within the twist-three framework, we can also derive the hard factors at leading order by following a similar analysis as in Ref. [26]. The method for calculating the single transversespin asymmetry for hard scattering processes in the twist-three approach has been developed in Refs. [43][44][45][51][52][53][54][55][56][57][58][59][60][61][62][63][64][65][66][67]. The collinear expansion is the central step to obtain the final results. We perform our calculations in a covariant gauge. The additional gluon from the polarized hadron is associated with a gauge potential A µ , and one of the leading contributions comes from its component A + . Thus, the gluon will carry longitudinal polarization. The gluon's momentum is dominated by x g P A +p g⊥ , where x g is the longitudinal momentum fraction with respect to the polarized proton. The contribution to the single-transverse-spin asymmetry arises from terms linear in p g⊥ in the expansion of the partonic scattering amplitudes. When combined with A + , these linear terms will yield ∂ ⊥ A + , a part of the gauge field strength tensor F ⊥+ in Eq. (9). Since p g⊥ = p ′ 2⊥ −p 2⊥ , the p g⊥ expansion of the scattering amplitudes can be performed in terms of the transverse momenta p 2⊥ and p ′ 2⊥ , which we can parametrize in the following way, The leading order diagrams shown in Fig. 2 can be calculated following the general procedure discussed above. The method is similar to the analysis of Drell-Yan lepton pair production in Ref. [65]. For example, to perform the Fourier transform of the SSA from transverse momentum to the impact parameter b ⊥ -space, we take the total transverse momentum of the two jets q ⊥ . In the leading order diagrams as shown in Fig. 2, the total transverse momentum q ⊥ can be easily identified: q ⊥ = p ′ 2⊥ for the left diagrams and q ⊥ = p 2⊥ for the right diagrams. Because the phases are opposite to each other for the left and right diagram, their total contribution will lead to the expression ǫ αβ S α ⊥ q β ⊥ = ǫ αβ S α ⊥ p β g⊥ . Using the definition of the twist-three matrix element, we find that the SSA contribution in the impact parameter b ⊥ -space can be written as . The hard factors can be calculated accordingly.

B. Soft Gluon Radiation
At one-loop order, we have to consider real gluon radiation associated with the production of the dijets. When the radiated gluons are parallel to the incoming partons' momenta, their contributions can be factorized into the associated parton distribution functions (from the unpolarized nucleon) or the polarized quark Sivers function (from the polarized nucleon). The gluon radiation will generate finite transverse momentum. According to the analysis of Ref. [26], we can write down the spin-dependent differential cross section as where P (<) represents the collinear splitting kernel excluding the end point contribution.
For the twist-three function it is given by A similar (albeit slightly simpler) expression holds for P (<) a ′ →a ⊗ f a ′ (x 1 ). In Eq. (31), the first term in the bracket comes from the collinear gluon radiation associated with the polarized nucleon, whereas the second term is associated with the unpolarized nucleon. An explicit calculation of all the relevant diagrams was presented in Ref. [26] for one particular channel (qq ′ → qq ′ ), and factorization arguments were given for all other channels. Only the so-called soft-and hard-gluonic poles are considered in the SSA calculations. However, all other pole contributions and T F contributions can be analyzed as well and similar results are expected.
In the following, we will focus on soft gluon radiation. The factorization of these contributions is more involved for several reasons. First, they contain double logarithms. In terms of transverse momentum distributions, we will find terms of the form 1/q 2 ⊥ ln(P 2 T /q 2 ⊥ ). These double logarithmic terms come from gluon radiation associated with all external particles. The collinear factorization arguments in Ref. [26] do not apply to these soft gluon emissions. Second, we have to deal with the soft gluon radiation associated with the final state jets. Using recent developments for the unpolarized case, we will be able to derive their contributions to the spin-dependent cross sections. We will first discuss several general features of twist-three calculations of the soft gluon radiation contributions to the SSA, and then we will apply these to the different partonic channels.

Generic Features of Twist-three Calculations
In Fig. 3, we show the generic diagrams that need to be calculated to obtain the soft gluon real radiation contributions. We follow the same strategy as in Ref. [17] to evaluate these diagrams. The radiated gluon carries transverse momentum k g⊥ which will contribute to the total transverse momentum of the two jets. The spin-dependent differential cross section can for a given partonic channel be schematically written as for the one-loop soft gluon radiation, where x 1,2 are defined as for leading order kinematics. This is because the soft gluons do not modify the longitudinal momentum fractions of the incoming partons. The partonic cross section H β twist-3 depends on the total transverse momentum q ⊥ , the hard momentum scale represented by P T (or in general,ŝ,t andû) and the jet radius R. Similar to the unpolarized case discussed in the last section, we have to exclude the gluon emission contributions belonging to the final state jets. Therefore, the partonic cross sections will depend on the jet size R.
From the diagrams in Fig. 3, we find that at this order the total transverse momentum of the two jets is equal and opposite to the transverse momentum of the radiated gluon: q ⊥ = − k g⊥ . Therefore, finite q ⊥ also implies finite k g⊥ . The main objective of the following calculations is to obtain H β twist-3 in the twist-three framework. Again, as briefly discussed above, we need to perform the collinear expansion for the incoming quark lines associated with the polarized nucleon. In the twist expansion, we take the limit of k g⊥ ≫ p 2⊥ ∼ p ′ 2⊥ ∼ p g⊥ . Meanwhile, we are also working in the correlation limit of q ⊥ ∼ k g⊥ ≪ P T . Therefore, the dominant contribution to the SSA comes from the expansion in powers of p g⊥ /k g⊥ . Any terms of the form p g⊥ /P T will be power suppressed in the correlation limit of q ⊥ ≪ P T .
To obtain the SSA for this process, the longitudinal gluon p g from the polarized nucleon needs to couple to the partonic scattering part to generate the necessary phase for a non-zero single spin asymmetry. Because we work at leading power in the limit of q ⊥ ≪ P T , we can classify the gluon attachments into two types. First, the p g gluon attaches to one of the initial/final state partons which does not radiate the soft gluon k g . Second, the p g gluon attachment and the soft gluon k g radiation happen on the same initial/final state parton. Here we discuss both cases, where the first type is easier to calculate, whereas the second type is somewhat more involved.
We first study the first type of diagrams. Particular examples are shown in Fig. 4, where the pole contribution comes from the gluon attachment to the final state parton k 1 and the soft gluon is radiated off the lines with momentum k 2 . Before an explicit evaluation, we would like to point out a number of important features which will help us to simplify the calculation. We will focus on the leading power contribution in the limit of small q ⊥ /P T . Therefore, the twist-three contributions only come from the p i⊥ -expansion associated with the radiated gluon line k g . For example, the p i⊥ dependence of the internal propagator (represented by the circle in Fig. 4) will lead to a power suppressed contribution in the limit of q ⊥ ≪ P T . Therefore, we only need to consider the p i⊥ -expansion in the propagators indicated by the red lines in Fig. 4. In addition, the lines cut by red bars are the places where we pick up the pole contributions. The left and right diagrams give opposite contributions from these two poles, because they are on opposite sides of the cut-line.
Because k 1 and k 2 are final sate observed momenta, it is convenient to keep them fixed in

FIG. 4:
Soft-gluonic pole contribution associated with the final state particle k 1 and the gluon radiation from the final state particle k 2 . The dashed line in the middle indicates the final state particles which are on mass-shell. The left diagrams represent the contribution when the gluon is attached to the left side of the cut-line, whereas the right diagrams correspond to the attachment to the right side of the cut-line.
the p i⊥ -expansion. As a consequence, the momentum flow will go through the radiated gluon momentum. For convenience, we define k g as the radiated gluon momentum with p i⊥ = 0, i.e., there is no p i⊥ dependence in k g . We label k gL as the momentum of the radiated gluon for the left diagram in Fig. 4, and k gR for the right diagram, respectively. Due to the fact that the momentum flows are different for these two diagrams, k gL and k gR will be different as well. Each of them is constrained by the on-shell condition for the radiated gluon. For example, we know that k gL⊥ = k g⊥ +p ′ 2⊥ . This gives the following momentum decomposition for k gL : We find that k 2 gL = 0 which satisfies the on-shell condition up to the linear term of p i⊥ . In the above expansion and the following calculations, we neglect all higher order terms of p i⊥ beyond the linear terms. Similarly, we find for the right diagram of Fig. 4, Once the kinematics are determined, we can proceed to calculate the soft gluon radiation contributions. This is similar to the unpolarized case. We multiply the Eikonal amplitudes of the diagrams shown in Fig. 4 and perform the collinear expansion of p i⊥ . For example, for the upper two diagrams, we have Right : where Γ µν (k g ) was defined in Eq. (14). We stress that k gL and k gR depend on p i⊥ , which implies that Γ µν will as well. Because the left and right diagrams give contributions with opposite sign for the phase, which is necessary to generate the SSA for this process, we will add their p i⊥ expansions with different signs. In the end, the total contribution from these two diagrams leads to the following expression: where we can simplify the above result as Fig. 4 Applying the twist-three procedure, the above leads to the following contribution to H β twist-3 , To determine the leading power contribution from soft gluon radiation, we need to integrate out the phase space of the radiated gluon. We will come back to this point after completing the analysis of all relevant diagrams. Now we turn to the lower two diagrams of Fig. 4. Here, k gL and k gR are the same as in the upper two diagrams. However, the p i⊥ -expansion comes from different propagators, Right : We also find that the final result is a little more complicated than for the upper two diagrams. After some algebra, we find that it can be written as The above derivations can be generalized to all other diagrams of the first type.
For the second type of diagrams, the longitudinal gluon from the polarized proton can attach to both the final state jet and the radiated gluon. Therefore, we have both soft-gluonic pole and hard-gluonic pole contributions. One particular example is shown in Fig. 5. Here, the gluonic pole contributions come from the final state particle k 1 which also radiated a soft gluon to generate the leading power contribution in the correlation limit of small q ⊥ /P T . The first diagram corresponds to the soft-gluonic pole and the rest to the hard-gluonic pole. The soft gluon pole leads to a delta function δ(x 2 − x ′ 2 ). In general, the hard gluon pole leads to a different delta function. However, in the correlation limit, the hard pole reduces to the same delta function. For example, the hard pole (labeled by the red bar in Fig. 5) leads to the following kinematics: In the correlation limit 1 , i.e., k µ g ≪ k µ 1 , the above reduces to x ′ 2 = x 2 . Therefore, the softgluonic and hard-gluonic pole contributions come from the same kinematics and there will be cancelations among them. These cancelations are very similar to those occurring for the SSA in the Drell-Yan process demonstrated in Ref. [52]. In particular, because of the color factor if abc T c = [T a , T b ], we can decompose the last diagram into the other two diagrams associated with the hard-gluonic pole. The combination exactly cancels out the soft-gluonic pole contribution from the first diagram. We have explicitly checked this cancelation for all these diagrams. To carry out the calculation, we have to follow the transverse momentum flow, and perform the p i⊥ expansion. The method is the same as that used to calculate the first type of diagrams: the kinematics of k g and p g will be determined from the on-shell conditions for k g and the pole contributions. More importantly, similar to the first type of diagrams, we only need to take into account the p i⊥ expansion from the denominators of the relevant propagators, since the expansions of the numerators are power suppressed. Therefore, the calculations are the same for all partons in the initial/final state, regardless of whether they are quarks or gluons. Both have the same denominators.
In the end, the total contribution will be a combination of the last two diagrams with the color factor of the third one. In summary, we can represent all four diagrams as the one on the left side. This can be repeated for diagrams associated with the initial state interaction with p 1 and the final state interaction with k 2 . We show the relevant diagrams in Figs. 6 and 7. These "effective" diagrams will be among the important ingredients for the final results. It is convenient to add the contributions from the first type and second type of diagrams. As an example, in Fig. 8, we show those diagrams for the initial state interaction contributions. These diagrams show that we can add soft gluon radiation on top of the initial state interaction diagram. The core part is the same for these three, in particular, the associated color factors. We will work out the color factors for the different channels later on. Here, we focus on the kinematics of the soft gluon radiation contribution and in particular on the leading power contributions.
As shown in Fig. 4, the contributions to the SSA come from the interference between the diagrams of Fig. 8 and the diagrams of Fig. 1. There will be diagrams where the longitudinal gluon is attached to the left side of the cut-line and to the right side of the cut-line. For convenience, we label these soft gluon radiation diagrams by their association with the external momenta. For example, we will label the first diagram of Fig. 8 by p µ 1 , the second diagram by k µ 1 and the third by k µ 2 . We label the diagrams in Fig. 1 similarly. The interference between the second diagram of Fig. 8  calculated above and the result is Similarly, we find the following result for the interference between the second diagram of Fig. 8 and the third diagram of Fig. 1: The calculation for the interference between the first diagram of Fig. 8 and the diagrams in Fig. 1 is much more involved. However, after a lengthy derivation, we find the results are very similar to the above two, and they can be all summarized as follows: Interestingly, we note that there is a one-to-one correspondence between the above results and those for the gluon radiation contributions for the unpolarized case in Sec. II. This is a very important feature to obtain the final factorization result for the SSA. The above analysis can be extended to diagrams with final state interaction contributions with k 1 and k 2 . For completeness, we show the diagrams in Fig. 9 and 10. The final results are the same as those given above for the initial state interaction contributions.
As mentioned above, in order to derive the leading power contribution to the SSA for this process from the above terms, we need to integrate over the phase space of the radiated gluon, where we keep the transverse momentum q ⊥ = − k g⊥ . Let us first work out the simple example of the p µ 1 p ν 1 term. This term is similar to that for Drell-Yan lepton pair production calculated in Ref. [52]. The phase space integral takes the following form where ξ = k g · p 2 /p 1 · p 2 and the lower limit of the ξ-integral comes from the kinematic limit ξ 0 = k 2 g⊥ /Q 2 . Working out the integral, we arrive at which is consistent with the double logarithmic behavior for Drell-Yan lepton pair production calculated in Ref. [52].
On the other hand, we can also carry out the above integral using integration by parts as follows: Of course, this gives the same result as above. However, this provides a convenient way to derive other terms as well. The rule is that the derivative only acts on the 1/k 2 g⊥ , not on the logarithmic terms. Taking the example of S g (k 1 , p 2 ) and S g (k 1 , p 1 ), we can follow the strategy and construct the following two terms: In the above two equations, R β does not contain a ln(1/k 2 g⊥ ) term. Therefore, we can perform the integration by parts directly and obtain the final result For T β , we notice that S g (k 1 , p 2 ) + S g (k 1 , p 1 ) = S g (p 1 , p 2 ) + 4 k 1⊥ · k g⊥ /(k 2 g⊥ k 1 · k g ), where the first term has been calculated above and the second term does not have a term ln(1/k g⊥ ). Applying this, we arrive at the following result for T β : Note that in practice we do the algebra and phase space integration in dimensional regularization in d = 4 − 2ǫ dimensions. For simplicity, we are not displaying here terms of O(ǫ). From the results for T β and R β , we are able to derive the corresponding results for the terms associated with S g (k 1 , p 1 ) and S g (k 1 , p 2 ). Following the same technique, we can also derive the results for S g (k 2 , p 2 ) and S g (k 2 , p 1 ). For S g (k 1 , k 2 ), since it does not contain a ln(1/k 2 g⊥ ) term, we can directly carry out the derivative with respect to k β g⊥ . Finally, we summarize all results here: Again, we have neglected the terms of O(ǫ) for simplicity. It is interesting to note that all of the above terms contribute to the leading double logarithms, whereas only k µ 1 k ν 1 and k µ 2 k ν 2 contribute to the jet related logarithms. Therefore, we need to consider all of them to derive the leading double logarithmic contributions. Next, we need to combine the above results with the associated color factors for the different channels in order to obtain the contributions to the SSA.
Let us first derive the SSA for the simplest channel, qq ′ → qq ′ , the quark-quark scattering with different flavors. This channel only has a t-channel diagram. The leading order results have been calculated in Ref. [26]. The hard factor is given by where the color factors for the initial and final state interactions are: For the initial state interaction contributions, we calculate the interference between the diagrams in Fig. 1 and Fig. 8. We obtain the following associated color factors: In order to obtain the final results, we multiply the leading power contributions of Eqs. (60)-(65) with the associated color factors. Adding these results, we obtain the leading contribution from soft gluon radiation which can be written as Here we only kept the terms relevant at LLA ′ , and we have h . We will come back to the remaining terms in Sec. IV when we discuss factorization breaking effects. The minus sign in the above equation is due to the fact fact that C I was computed on the basis of the quark Sivers function for the SIDIS process, which has an opposite sign compared to the normalization of the twist-three matrix element T F . It appears that the terms in Eq. (73) do have a clear factorization structure that includes the leading double logarithmic term and the terms associated with the final state jets. The latter are represented by logarithmic terms of the jet radii R 1,2 .
We can extend the above calculations to the final state interaction contributions associated with k 1 and k 2 . We summarize the relevant color factors for the different terms in Table I. The total contribution can be written as where H Sivers qq ′ →qq ′ was defined in Eq. (66). Clearly, the first term contributes to the leading double logarithms of the SSA for this process. In addition, the divergence associated with the two final state quark jets has the desired structure. Combining the above result with the collinear gluon radiation contributions from incoming partons, see Eq. (31), we obtain the spin-dependent differential cross section for the qq ′ → qq ′ channel in the correlation limit of q ⊥ ≪ P T . When taking the Fourier transform to b ⊥ -space and adding the virtual and jet contributions to cancel the divergences, we expect to obtain the following one-loop result for W T β at LLA ′ accuracy: Here we have included the subtraction of the collinear divergences, and P T qg→qg and P a ′ →q are the complete splitting kernels for the associated twist-three and leading-twist parton distributions, respectively. To derive the above results, we have assumed that the virtual contributions cancel the soft divergences of the real gluon radiation. It is important to show that this is indeed the case. However, this is beyond the scope of this work and we plan to come back to this issue in a later publication.
We summarize three important aspects of the above result at this order. First, the collinear splitting is associated with the twist-three and twist-two parton distribution functions. This is the essence of the factorization of collinear gluon radiations as demonstrated in Ref. [26]. Second, the result for the leading double logarithms is consistent with the collinear and soft factorization at LLA ′ . Each of the incoming quark lines contributes half of this double logarithmic term. This is an important feature for the Sudakov resummation. Third, the logarithms associated with the jets are also factorized in terms of the individual jets and we obtain the expected color charges of the final state jets. These results are consistent with the factorization argument.
In the following two subsections, we will extend the above analysis to two other important channels, gq → gq andqq → gg. The former is the dominant channel for the SSA in dijet production for the typical kinematics at RHIC.

qg → qg
The leading order derivation for the channel qg → qg was carried out in Ref. [26]. In order to simplify the analysis for the soft gluon radiation contribution, we follow Ref. [17] and decompose the fundamental partonic scattering amplitude as with two different color structures at the Born level. Here, a and b are the color indices for the incoming and outgoing gluons, and i and j for the incoming and outgoing quarks. The amplitudes A 1 and A 2 depend on the momenta of the two incoming particles, p 1 and p 2 , and on the momenta of the two outgoing particles k 1 and k 2 for the quarks and gluons, respectively. The single spin asymmetry can be formulated starting from the decomposed amplitude above. For example, the initial state interaction contribution can be derived from the following expression: The color indices i and i ′ are coupled to the adjoint representation of the twist-three Qiu-Sterman matrix element. Therefore, we can rewrite u iūi ′ as T c ii ′ , and the above result leads to where N color = (N 2 c − 1)N c C F . Similarly, for the final state interaction with the gluon line, we have For the final state interaction contribution from the final state quark line, we have Noticing that the initial state interaction carries a minus sign, the total contribution will be From our parameterization of the Born amplitude, we have Substituting the above expressions into Eq. (82), we reproduce the result for the leading order H Sivers qg→qg in Ref. [26]. Now, let us turn to the soft gluon radiation contributions. For the initial state interactions, we have soft gluon radiation contributions from the incoming gluon line, the outgoing gluon line and quark line. The relevant diagrams will follow those in Fig. 8. In this case, the upper two lines are gluons. The associated amplitudes for these three diagrams are given by Here c is the color index of the radiated gluon and d corresponds to the longitudinal gluon from the polarized nucleon attached to the partonic scattering part. The contributions to the SSA come from the interference of the above amplitudes and those in Fig. 1, which we list here for the qg → qg channel, Similar to the previous case, the following interference terms are simple, which will be multiplied by the leading order initial state interaction contribution of Eqs. (60)- (62). The other interference terms are a little bit more complicated. We have By combining them with the associated kinematic contributions in Eqs. (60)-(65) and adding them, we obtain the leading logarithmic result for the SSA from the initial state interaction for this channel. In particular, the above three terms add up to which exactly cancels the leading logarithmic contribution from the p µ 1 p ν 1 term which also has the color factor C A . We thus obtain the final result for the initial state interaction contribution: For the final state interaction associated with the final state gluon jet, we have the following amplitudes for the associated diagrams in Fig. 9: Again, the interference contributions from p µ 1 p ν 1 , k µ 1 k ν 1 and k µ 2 k ν 2 have the same structure as those for the initial state interaction diagrams in Eq. (91). The remaining contributions can be written as follows: Adding up the three contributions, we have The total contribution from the final state interaction with the gluon jet is given by For the final state interaction with the quark line, we find the following amplitudes for the associated diagrams of Fig. 10: and the interference terms are The above three terms lead to The final result from the final state interaction with the quark line is Adding up all initial and final state interaction contributions, we obtain the following result for the SSA at the leading logarithmic order: Including the collinear gluon radiation contributions from the incoming partons as in Eq. (31), we obtain the spin-dependent differential cross section for the gq → gq channel in the correlation limit of q ⊥ ≪ P T which is given by After taking the Fourier transform to b ⊥ -space, we expect to have the following one-loop result for W T β : The computation for the qq → gg channel is very similar to the case discussed above. Here, we can parametrize the leading Born amplitude as where, of course, the amplitudes will be different from those for the qg → qg channel.
In particular, we have ( for the current case.
The single spin asymmetry comes from the initial state interaction with the antiquark line and the two final state gluon lines. For the initial state interaction contribution, we have where N (qq) color = N 2 c C F . For the final state interaction contribution with the gluon line of k 1 , we obtain C (qq) Using the the symmetry of the channel, the final state interaction contribution from the gluon line k 2 can be obtained from the above result as The total contribution at the leading order is given by Substituting the expressions of A 1 and A 2 into the above equation, we reproduce the result for the leading order H Sivers qq→gg in Ref. [26]. Next, we can work out the soft gluon radiation contribution as well. For the initial state interaction contributions (C I ), we consider the soft gluon radiation from the incoming antiquark line, and two the outgoing gluon lines. The relevant diagrams again correspond to those shown in Fig. 8. The associated amplitudes for these three diagrams are where c is again the color index of the radiated gluon and d corresponds to the longitudinal gluon from the polarized nucleon. The single spin asymmetry contributions come from the interference of the above amplitudes and those of Fig. 1. For qq → gg channel, we have In our previous notation, the resulting interference terms are The other interference terms can be evaluated as The above three terms add up to which exactly cancels the leading log contribution from k µ 1 k ν 1 and k µ 2 k ν 2 with color factor C A . Therefore, we obtain the following final result for the initial state interaction contribution: For the final state interaction associated with the final state gluon jet k 1 , we find the following amplitudes, see Fig. 9, Again, the interference contributions from p µ 1 p ν 1 , k µ 1 k ν 1 and k µ 2 k ν 2 have the same structure as those for the initial state interaction diagrams. The remaining contributions can be written as Adding the three terms above, we find The total contribution from the final state interaction with the gluon jet leads to Again, from the symmetry under interchange of the final state particles, we can obtain the soft gluon radiation contribution for the final state interaction with the gluon line k 2 . Adding all these contributions, we have the following result for the SSA at the leading logarithmic order, Adding the collinear gluon radiation contribution, we find the following result: for the qq → gg channel. Taking again the Fourier transform to b ⊥ -space should lead to the following one-loop result: C. Factorization and Resummation at LLA ′ Our above results for the soft gluon radiation contribution at one-loop order are very instructive for developing a factorization argument according to which one can perform the resummation of logarithms in q ⊥ /P T . Three important types of contributions are incorporated and explicitly presented in our calculations: (1) collinear gluon radiation from the incoming hadrons, (2) collinear gluon radiation from the final state jets, and (3) the leading double logarithms from soft gluon radiation.
We argue that these contributions to the LLA ′ will continue to factorize at higher orders. As illustrated in Fig. 11, we can summarize the above calculations of the one-loop soft gluon radiations in a factorization structure. First, the power counting analysis simplifies the kinematics. As shown in this figure, the initial/final state interactions needed to generate a phase for the SSA in this process can be factorized into a hard partonic scattering amplitude, in analogy to the leading order Born amplitude for the spin-averaged case. Additional soft gluon radiation only appears in association with the external lines. This not only simplifies the detailed derivations at this order, but also shows that the soft gluon radiation associated with the final state jets can be generalized to all orders.
In our example, we have chosen a physical polarization for the radiated gluon, such that the jet contribution comes only from the squared diagrams where the radiated gluon is attached to one particular jet, e.g. either k 1 or k 2 . These emissions always result in a leading contribution proportional to 1/q 2 ⊥ ln(1/R 2 ). An evolution equation can be derived to resum higher order emissions and the final result can be written in terms of the Sudakov resummation form factor of Eq. (8).
Similar arguments can be made for the collinear gluon radiation associated with the incoming hadrons. These collinear gluon also contribute terms of order 1/q 2 ⊥ , which will be multiplied by the splitting of the associated parton distribution and/or twist-three Qiu-Sterman matrix element. The all order resummation can be carried out by solving the relevant DGLAP equations. This can be achieved by evaluating the distributions at the scale µ b = b 0 /b ⊥ and by including the associated anomalous dimensions in the Sudakov resummation form factor.
The leading double logarithms from soft gluon radiation are associated with kinematics overlapping with the collinear gluon radiation from the incoming partons. The resummation of these double logarithms can be carried out by solving the associated Collins-Soper evolution equations for the TMD parton distributions. The double logarithms from two TMD parton distributions give the leading double logarithms for the final differential cross section. This argument has been verified explicitly in the above derivation. We expect that this can be generalized to higher orders as well and an all order resummation can be obtained in the form of the Sudakov form factor in Eq. (8).

IV. FACTORIZATION BREAKING EFFECTS AT NLL
To exhibit factorization breaking effects beyond the LLA, we consider the channel qq ′ → qq ′ as an example. From the derivation in the last section, we obtain the soft gluon radiation contribution to the transverse spin dependent differential cross section (see Eqs. (60)-(65), (73), (74)) where h (0) Comparing to Eq. (66) of Ref. [17], we observe that this term is different from the analogous term in the unpolarized case, implying that it cannot be factorized into a spin independent soft factor for the qq ′ → qq ′ process. Even if we consider only initial state soft gluons for the transverse spin dependent differential cross section, we get a result different from the unpolarized one. Here we have where The same findings occur for the other two channels studied in the previous section, qg → qg and qq → gg. Let us recall the factorization argument for soft gluon radiation in dijet production. Following the analysis of the generic soft gluon radiation in this process, see, e.g., in Refs. [19][20][21], one can derive the associated soft factor for the TMD factorization in matrix form of the relevant color spaces [16,17]. As shown in Fig. 12, the color configuration in the 2 → 2 partonic process can be written as ij → kl, where ij and kl represent the color indices for the incoming and outgoing partons, respectively. The independent color space tensors carry The contributions arise as the interference between the amplitude for ij(a) → kl with an additional gluon attachment from the polarized nucleon and the amplitude without the gluon attachment i ′ j → kl, where a represents the color index of the attached gluon. The final result is derived by contracting color indices i and i ′ with T a ii ′ . Therefore, the color flow is totally different from that in Fig. 12 for the unpolarized case.
these indices. The relevant amplitudes for soft gluon radiation and their complex conjugates are derived based on these configurations as well. In TMD factorization, the associated soft factors are formulated in terms of a matrix form on the basis of these color configurations as well, where the external lines will be replaced by the Wilson lines.
However, for the SSA in the dijet process, an additional gluon attachment from the polarized nucleon is needed to generate the phase required for the Sivers effect. Because of this, the color flow of the partonic process is totally different from that in the unpolarized case. As shown in Fig. 13, the single spin asymmetry comes from the interference between the amplitude with gluon attachment and that without gluon attachment. The attaching gluon carries color (with adjoint representation index a in this diagram). The color flow for the amplitude shown in this plot can be written as ij(a) → kl for the left side and i ′ j → kl for the right side, respectively. The final result for the SSA is derived by contracting i and i ′ via the matrix T a ii ′ representing the attaching gluon a. Because of this additional color flow caused by the gluon attachment, the color space configurations differ from those in the unpolarized case. As a consequence, the soft factor changes, and the result for the spinaveraged case derived in Ref. [16,17] does not apply here. This is indeed precisely what our calculations demonstrate.
As just described, the difficulty to apply the standard soft factor factorization for the SSA comes from the fact that we need to generate a phase in order to obtain a nonzero SSA for this process. The phase comes from the imaginary part of the propagator (pole contribution), which may carry a different sign compared to the real part. For unpolarized scattering at this order, in a similar diagram like Fig. 13, the attached gluon can be incorporated into a relevant gauge link in the unpolarized quark distribution. This is only possible when there is no pole contribution from the propagators associated with the attached gluons. After this, the color index for the quark line entering the hard part from the left side will be the same as the on right side, i.e., i → i ′ . As a consequence, we can derive the associated soft factor in TMD factorization as indicated in Fig. 12. Of course, as mentioned in the Introduction, at higher order O(α 3 s ) and beyond, the factorization for the unpolarized cross section also breaks down. The reason is the same: one can have two poles in the propagators associated with two gluon attachments, which can not be factorized into the relevant gauge link in the TMD quark distributions; see, examples given in Refs. [27,30].
To summarize, the soft gluon radiation contribution to the SSA in the dijet production process in pp collisions cannot be factorized into a spin independent soft factor because of the pole contributions, thus breaking standard TMD factorization. Since the soft factor carries a next-to-leading logarithmic contribution, the factorization for the SSA breaks down at NLL as explicitly shown by the above example of qq ′ → qq ′ channel. At the leading logarithmic level, we can argue that the relevant soft gluon radiation belongs to the parton distributions and the final state jets. The associated large logarithms can be resummed through the evolution of these parton distributions and jet functions, and factorization is only broken beyond leading logarithmic accuracy. Therefore, the SSA for dijet production in polarized pp collisions may provide a unique opportunity to explore factorization breaking effects. We will study the associated phenomenology in the following section.
We note that it is conceivable that the soft gluon radiation in Fig. 13 may be factorized into a "non-traditional" soft factor that is unique to the SSA. This might be achieved by setting up different color tensors on the left side and the right side of the diagram. On the right side of the diagram in Fig. 13 we have a standard 2 → 2 partonic channel with its simple and standard color flow. On the left side, we have a 3 → 2 partonic channel, for which recent developments on the soft gluon evolution in 2 → 3 processes [69,70] may become useful. The color-contraction of the two sides is expected to give rise to a non-square matrix problem. It appears likely that the structure of one-loop soft gluon radiation in the singletransverse spin cross section may be understood in this way; whether one can generalize this to higher orders in TMD factorization would be an open question. It is worthwhile to further pursue a study along this direction, which, however, is beyond the scope of the current paper. We hope to address this in a future publication.

V. PHENOMENOLOGICAL STUDY AND COMPARISON TO STAR DATA
Within the improved leading logarithmic approximation, we can write the differential cross sections for the spin averaged and spin dependent cases in the correlation limit q ⊥ ≪ P T as follows: Following our arguments on the improved leading logarithmic factorization, the resummation formulas for the above W uu and W T β can be written as where was defined in Eq. (8). For the unpolarized case, we could in principle also use a more advanced resummed result as derived in Ref. [17]. For the following studies, we have decided to use the same resummation accuracy for both the polarized and unpolarized cross sections. We have checked numerically that this does not introduce any sizable effects for the unpolarized case as far as the spin asymmetries are concerned.
In order to compare to the experimental data, we have to make further approximations. This is because the phenomenology of the Sivers function (or the Qiu-Sterman matrix element) in the spin-dependent cross section is presently not sophisticated enough to warrant the inclusion of detailed evolution effects in our calculations. Therefore, we will apply estimates of the quark Sivers functions from known experiments without considering their scale dependence. We approximate the parton distributions in the above equations at a fixed lower scale, e.g., µ b ≈ µ 0 = 2GeV, where the quark Sivers functions are extracted from SIDIS. We comment on the uncertainties of these extractions below.
With these modifications, the resummation formulas can be written as where µ 0 will be varied around 2 GeV in our final results. In addition, these modifications also allow us to simplify the numerical calculations. Within the above approximation, the q ⊥ -dependence only comes from the Sudakov form factor (with the additional b β ⊥ for the Sivers effect) which also depends on Q 2 .

A. Preliminary results
According to the above results, the Sivers asymmetries depend on three ingredients: (1) the Sivers functions; (2) the associated hard factors for the different partonic channels; (3) the q ⊥ dependence from the Sudakov form factor and the related non-perturbative TMDs. In the following, we will briefly go through these three ingredients, and we then address the comparison to the STAR data [2,46].
FIG. 14: The Sivers functions used in our calculations. We use the extractions from Ref. [66] which are obtained from Sivers SSAs in semi-inclusive hadron production in deep-inelastic scattering.

Sivers Functions
The quark Sivers functions have been mainly determined from single transverse spin asymmetries in semi-inclusive hadron production in the deep-inelastic scattering process. The latest global analyses can be found in Refs. [71,72] (see also references therein). The associated Qiu-Sterman matrix elements T F (x, x) have also been extracted from the single spin asymmetry A N for single inclusive hadron production in pp collisions [54]. However, the latter process contains the Collins twist-three fragmentation contributions as well [73][74][75][76][77]. Therefore, it may not be sufficient to constrain the Sivers functions from inclusive hadron A N in pp collisions alone. Recently a first attempt was made to perform a global analysis of SSA data from different processes including A N in pp collisions [71].
In the following calculations, we will use the Sivers functions constrained by SSA data in SIDIS. We would like to add a few comments concerning the precision of these extractions. First, the current experimental data on SSAs in SIDIS come from the relatively low Q 2 region, where TMD factorization may not be as rigorous as compared to high Q 2 . This situation will of course be improved by the future Electron-Ion Collider. Until then, we have to keep in mind the systematic uncertainties of the Sivers functions extracted from the existing SIDIS data. A crosscheck with other processes, such as Drell-Yan lepton pair production and W /Z-boson production in pp collisions will provide crucial information on the Sivers functions. Second, the quark Sivers functions determined from SIDIS have a particular feature -the up quark and down quark distributions have opposite signs. As a result, one often finds a strong cancelation between these two quark Sivers functions when they are used in a physical process, resulting in significant uncertainties of phenomenological extractions of the functions. All existing fits report at least 10% uncertainty of the valence quark Sivers functions. In addition, the sea quark Sivers functions have even larger uncertainties. We hope that the future EIC will help to better constrain both valence and sea quark Sivers functions.
For our numerical results we will use the quark Sivers function from Ref. [66] as a baseline, keeping in mind their significant uncertainty just described. The dijet asymmetries studied in this paper are precisely a case where the u and d quark Sivers functions are added (at least for the dominant qg channel) and cancel to a significant extent. This also increases the uncertainty of any predictions that can be made. Given the uncertainty of the u + d combination, and to obtain a conservative order of magnitude estimate, we assume the total u + d distribution to be ±0.2u, allowing the distribution to have either sign. We also neglect the sea quark Sivers function contributions in the numeric estimate. Because of the qg → qg channel dominance, the individual difference between the u and d quark Sivers in our parameterization will not affect much the final results.
More recently, the STAR collaboration has developed a novel approach to study the SSA in dijet production at RHIC using quark flavor tagging by measuring the charge of jets [46]. By tagging the total charge of the final state jet, one can separate u-quark jets from d-quark jets which potentially avoids the cancelation between the two distributions to some degree. The SSAs for the flavor tagged dijet production will be directly related to either the u-quark Sivers function or the d-quark Sivers function, depending on the charge of the triggered jet. We will comment on the comparison to these exciting new data at the end of this section. See also Ref. [78] for the original proposal of the jet charge by Feynman and Field, Refs. [79,80] for recent theoretical studies and Refs. [81,82] for experimental measurements by ATLAS and CMS.

Hard Factors
We expect that the quark-gluon channel will make the dominant contribution to the SSA in dijet production, especially at forward rapidity where one probes the valence region of the Sivers functions, which in fact is primarily constrained by the SIDIS experiments. It is instructive to examine the hard factor for this channel as an example to study the kinematic dependence of the SSA of our process.
In Fig. 15, we show the ratio of the spin-dependent hard factor and the spin-averaged hard factor for the qg → qg channel for typical kinematics at RHIC. The leading jet transverse momentum is chosen as P T = 6 GeV and the rapidity interval we consider is −1 < y 1 < 2. From the figure, we can clearly see that the ratio increases with rapidity. In particular, the asymmetry will be larger in the forward region.
Another important feature of the hard factor is that it is positive for all kinematics. This means that the asymmetry is dominated by final state interaction effects. This is different from the inclusive hadron A N , which is dominated by initial state interaction effects. The reason is that for dijet production, the final state interaction with the gluon line in the qg → qg channel cancels the initial state interaction with the initial gluon line. On the other hand, for single inclusive hadron production, the final state interaction with the gluon line does not contribute to the quark fragmentation part in the qg → qg channel which is the dominant source for hadron production in pp collisions. Therefore, the dijet SSA will have an opposite sign compared to the single inclusive hadron SSA. This is a very interesting feature and will have significant phenomenological consequences for SSAs at RHIC.

Sudakov Effects
It has been known for some time that the Sudakov effects will suppress the single spin asymmetries [83]. This suppression factor was also included when the dijet SSA was proposed in Ref. [25]. Here, we would like to follow up on this issue and discuss the effect on the dijet asymmetries using the updated results for the Sudakov form factor and non-perturbative TMDs. We will continue to focus on the qg process.
As mentioned at the beginning of this section, the q ⊥ -dependence is contained in the Sudakov form factor and the associated non-perturbative TMDs. For the discussion here and the numerical results, we separate the q ⊥ -dependence of the spin-dependent and spinaveraged differential cross sections as where J 0,1 are Bessel functions and where we have applied the b * -prescription in the above equation The perturbative part of the Sudakov form factor is the same for both cases S (qg) pert (Q 2 , b * ) = S T (qg) Sud (Q 2 , b * ) which was given in Eq. (8). The non-perturbative parts are parametrized as [84] S (qg) with g 1 = 0.212, g 2 = 0.84, g s = 0.062 and Q 2 0 = 2.4 GeV 2 . As an example, we plot in Fig. 16 the ratio R T /R as a function of q ⊥ for typical values of Q relevant for RHIC kinematics. Clearly, the asymmetry increases with q ⊥ . Different from previous examples (SIDIS or Drell-Yan lepton pair production), the curves in the plot do not reach the maximum of the asymmetry in the q ⊥ region relevant for TMD physics. The reason is that for the qg channel, the Sudakov form factor leads to a strong q ⊥ -broadening effect. In particular, this is due to the the double logarithms associated with the incoming gluon distribution which push the peak of the asymmetry to higher values of q ⊥ , beyond the TMD domain.  16: The ratio of R T (Q, q ⊥ ) and R(Q, q ⊥ ) as a function of q ⊥ for typical Q values to dijet production at RHIC.

B. Comparison to the STAR Data from 2007
As suggested in Ref. [25], the dijet spin asymmetry can be measured through the azimuthal angular distribution between the two jets. Because the Sivers asymmetry leads to a preferred direction of the total transverse momentum of the two final state jets, the angular distribution will be shifted toward that direction related to the Sivers asymmetry. The magnitude of the asymmetry will depend on the relative angle between the leading jet and the polarization vector S ⊥ . In particular, as mentioned in Ref. [25], the SSA for dijet production is at its maximum value when the jet direction is parallel or anti-parallel to the spin vector S ⊥ of the proton. However, the asymmetry will vanish if the leading jet is perpendicular to the spin S ⊥ . It can be shown that this introduces an additional factor of | cos(φ j − φ S )| for the dijet SSA. Therefore, when we compare to the STAR data, we need to include an average of this factor over the azimuthal angle between the leading jet and the polarization vector: 1 π π 0 dφ| cos(φ)| = 2 π . As mentioned above, for the dijet SSA we take the u+d Sivers distribution to be 20% of the extracted u-quark Sivers function. From the existing experimental data, we can determine neither the size nor the sign of the total contribution from the u and d-quark Sivers functions. Therefore, to estimate the total contributions to the dijet SSA, we will use both signs, in order to obtain an idea of the uncertainties associated with these determinations.
The jet kinematics of the data published by the STAR experiment in 2007 is P T > 4 GeV and rapidity in the range of −1 < y 1,2 < 2. In order to compare to the experimental data, we integrate out the leading jet momentum and the relative rapidity between the two jets, but we keep the total rapidity y = y 1 + y 2 .
Using the differential cross sections for the spin-averaged and spin-dependent cases in the MLLA given above, we obtain the following expression for the single transverse spin asymmetry which can be compared to the STAR measurement: A N (y) = 2 π acd b=u,d d 2 P T dy 1 x 1 f a (x 1 , µ 0 ) x 2 T F b (x 2 , µ 0 ) H Sivers ab→cd (P T , x 1 , x 2 ) w T (Q) abcd d 2 P T dy 1 x 1 f a (x 1 , µ 0 ) x 2 f b (x 2 , µ 0 ) H uu ab→cd (P T , x 1 , x 2 ) w(Q) , where Q 2 =ŝ = x 1 x 2 S pp and w(Q), w T (Q) are weights for the total transverse momentum The SSA in dijet production at RHIC as a function of the total rapidity y = y 1 + y 2 of the two jets for the kinematics of the STAR measurement in 2007: P T > 4 GeV and −1 < y 1,2 < 2.
The upper bound corresponds to T F u +T F d with +20% of the fitted value of T F u of Ref. [66], whereas the lower bound corresponds to −20%. integral, The upper limits of the above integrals correspond to the TMD region where q ⊥ ≪ P T . Notice that Q ≥ 2P T for all kinematics. With our assumptions on the u and d Sivers functions, we calculate the SSA in Eq. (158), and find that the asymmetry is of order 10 −4 for the entire rapidity range shown, see, Fig. 17. We note that the central values of u-quark and d-quark Sivers functions from the fit of Ref. [66] lead to smaller asymmetries shown in Fig. 17. Let us summarize the main differences with respect to the previous calculation in Ref. [31]: First, the Sivers functions determined from SIDIS are different. Second, we have included the relative suppression from Sudakov effects. If we integrate over the entire rapidity region, the asymmetry is about 1.7 × 10 −4 .
It is interesting to note that the STAR measurement in 2007 found that the SSA for dijet production is consistent with zero, and their uncertainties are of the order of 10 −2 . According to our calculation, a finite asymmetry could be observed if the uncertainty is reduced by more than one order of magnitude. Of course, this also depends on the size of the total up and down quark Sivers function. If they completely cancel, then the asymmetries will depend on the sea quark Sivers functions, which are known to be not as large as the valence ones.

C. The Flavor Tagged Dijet Asymmetry
In the last few years, the STAR collaboration has investigated a novel method to explore the SSA in dijet production. They considered the jet charge to tag the up or down quark flavor of the jet. An up (down) quark jet has positive (negative) jet charge, whereas a gluon jet leads to a neutral jet charge. From the preliminary analysis, the efficiency of this separation is reasonably good. A similar idea is also proposed in [85]. This also suggests further improvements by using the jet flavor information in the jet charge definition.
In order to compare to the experimental data, we take into account the u and d quark contributions separately, just for the qg → qg channel. For up quarks, we have A (up) N (y) = 2 π d 2 P T dy 1 x 1 f g (x 1 , µ 0 ) x 2 T F u (x 2 , µ 0 ) H Sivers gq→gq (P T , x 1 , x 2 ) w T (Q) d 2 P T dy 1 x 1 f g (x 1 , µ 0 ) x 2 f u (x 2 , µ 0 ) H uu gq→gq (P T , x 1 , x 2 ) + (x 1 ↔ x 2 ) w(Q) . (161) An analogous expression holds for the d-quark Sivers contribution. In Fig. 18, we plot the two asymmetries as functions of y = y 1 + y 2 . From this plot, we can see that the asymmetries are of the order of 10 −3 . If we integrate over the entire rapidity range, we obtain an asymmetry of 2.2 × 10 −3 and −3.7 × 10 −3 for the ug → ug and dg → dg channels, respectively. Our results compare reasonably well to the preliminary STAR results 2 .
Compared to the results shown in Fig. 17, we find that the asymmetries are much larger for the flavor tagged case. This is not only due to the fact that for the flavor tagged case no cancelation between the u and d-quark Sivers functions occurs, but also because the denominator of the unpolarized cross section only contains the qg → qg channel. By tagging the (quark) flavor of the final state jets, we exclude a major background contribution from gg → gg, which only leads to charge neutral jets in the final state.
The asymmetries shown in Fig. 18 assume 100% efficiency of the tagging in the jet sample. To compare to the STAR data properly, we need to consider a realistic tagging efficiency, which will suppress the asymmetries somewhat.

VI. SUMMARY AND DISCUSSION
We have investigated the single transverse spin asymmetry in dijet correlations in hadronic collisions. The total transverse momentum of the dijet in the final state is correlated with the incoming nucleon's polarization vector. We have focused on the SSA contribution from the quark Sivers function of the polarized nucleon where initial and final state interaction effects play an important role. A detailed analysis at one-loop order has been carried out for the contribution from soft gluon emissions in order to understand the factorization properties. It was found that the associated TMD factorization is valid at the level of leading double logarithms and single logarithms from the TMD quark distribution and those collinear to the jet. However, additional soft gluon contributions to the single logarithms do not show the same pattern as in the unpolarized case investigated in Ref. [17] and hence cannot be incorporated in the TMD factorization formula in terms of the same spin independent soft factor. This indicates that the standard TMD factorization is broken at the single logarithmic level for the SSA in dijet correlations in hadronic collisions. We believe that this issue will deserve further attention by investigating whether a consistent "non-traditional" soft factor for the single transverse spin asymmetry could be defined.
We have further presented phenomenological studies for the SSA in dijet production at RHIC based on the LLA ′ approximation, for which one improves the standard LLA resummation formula by "universal" subleading logarithmic terms associated with the initial partons and the final state jets. Using the quark Sivers functions constrained by SSAs in SIDIS, we have found that the leading double logarithmic resummation effects suppress the asymmetry for the kinematics relevant for the measurements by the STAR collaboration at RHIC, making them broadly consistent with experimental results. We also presented results for the flavor (charge) tagged dijet case, where the asymmetries are much larger than when the jet flavor is not tagged. A detailed comparison with the experimental data will be helpful to understand factorization breaking effects.
We finally note that in our analysis we have only considered the perturbative part of the factorization breaking effects. The non-perturbative TMD quark distribution could also contribute to such effects. Our numerical results assume that this part is the same as for SIDIS. To address this issue, more detailed comparisons to experimental measurements and further theoretical efforts are necessary. In any case, a more refined phenomenology will be warranted once there is a better general understanding of factorization breaking effects, along with more detailed experimental information. wherex = k g · p 1 /p 1 · p 2 . Clearly the collinear divergence associated with the jet is canceled between these two terms. Averaging over the azimuthal angle of the jet, we have where ξ 0,1 are the same as above. Carrying out the integral, we find Together with the result for S g (p 1 , k 1 ) we obtain the result for S g (P, k 1 ). A similar result can be obtained for S g (P, k 2 ), and we are able to derive the result for S g (k 1 , k 2 ) in Eq. (27).