The Sivers Asymmetry in Hadronic Dijet Production

We study the single spin asymmetry in the back-to-back dijet production in transversely polarized proton-proton collisions. Such an asymmetry is generated by the Sivers functions in the incoming polarized proton. We propose a QCD formalism in terms of the transverse momentum dependent parton distribution functions, which allow us to resum the large logarithms that arise in the perturbative calculations. We make predictions for the Sivers asymmetry of hadronic dijet production at the kinematic region that is relevant to the experiment at the Relativistic Heavy Ion Collider (RHIC). We further compute the spin asymmetries in the selected positive and negative jet charge bins, to separate the contributions from $u$- and $d$-quark Sivers functions. We find that both the sign and size of our numerical results are roughly consistent with the preliminary results from the STAR collaboration at the RHIC.


I. INTRODUCTION
Exploring transverse momentum dependent parton distribution functions (TMD PDFs) has become one of the major research topics in hadron physics in recent years [1]. TMD PDFs provide three-dimensional (3D) imaging of the nucleon in both the longitudinal and transverse momentum space, which is one of the scientific pillars at the future Electron-Ion Collider [2]. Such 3D imaging of the nucleon offers novel insights into the highly nontrivial non-perturbative QCD dynamics and correlations [3].
Sivers function is one of the most studied TMD PDFs in the community. It describes the distribution of unpolarized partons inside a transversely polarized nucleon, through a correlation between the transverse spin of the nucleon and the transverse momentum of the parton with respect to the nucleon's moving direction. The Sivers function was first introduced by Sivers in 1990s [4,5] to describe the large single transverse spin asymmetries observed in single inclusive particle production in hadronic collisions, see e.g. [6,7]. Since then, large single spin asymmetries have also been consistently observed in proton-proton collisions in high energy experiments at the Relativistic Heavy Ion Collider (RHIC) [8][9][10][11][12][13]. On the theoretical side, understanding the precise origin of such large spin asymmetries has triggered extensive research in the QCD community [14][15][16][17][18][19][20][21][22]. The difficulty in understanding such asymmetries for single hadron production (such as pions) in proton-proton collisions lies in the fact that they could receive contributions from many different correlations. Besides Sivers type correlations, whose collinear version is referred to as the Qiu-Sterman function [15,23] in the incoming nucleon, there could also be similar correlations in the hadronization process when the parton fragments into the hadrons [17,19,21,22,24]. See [25] for a recent development along this direction.
Simultaneously the Sivers asymmetry has also been studied in semi-inclusive deep inelastic scattering (SIDIS) by HERMES collaboration at DESY [26,27], COMPASS collaboration at CERN [28,29], and Jefferson Lab [30]. Because of the semi-inclusive nature of the process, one can isolate the contribution from the Sivers function via different azimuthal angular modulations [31]. One of the remarkable and unique properties of the Sivers functions is its non-universality nature. For example, based on parity and time-reversal invariance of QCD, one can show that quark Sivers functions in SIDIS are opposite to those in the Drell-Yan process [32][33][34]. Such a sign change has been studied and confirmed experimentally [35][36][37][38], though additional work remains to be done to quantify the change in more details [39].
Sivers effect has been continuously studied in proton-proton collisions at the RHIC. In order to eliminate the contributions from the spin correlations in the fragmentation process, the Sivers asymmetry for jet production processes has been explored in the experiment [13,40,41]. In particular, back-to-back dijet production in transversely polarized The rest of the paper is organized as follows. In Section II, we summarized the QCD formalism for dijet production in both unpolarized and polarized scatterings, and we provide a few remarks about our formalism. In Section III, we provide a procedure and demonstrate how to compute the process-dependent polarized hard functions in the color matrix form. In Section IV, we present the renormalization group evolution of all the relevant functions in our formalism, and we provide the final resummation formula. Section V is devoted to the phenomenological studies, where we make predictions for dijet Sivers asymmetry in the kinematic region relevant to the experiment at the RHIC. Since we are mainly interested in the Sivers asymmetry in the forward rapidity region where quark contributions dominate, we consider only the quark Sivers contribution and neglect the gluon Sivers contribution. We summarize our paper in Section VI.

II. QCD FORMALISM FOR DIJET PRODUCTION
In this paper, we study back-to-back dijet production in transversely polarized proton-proton collisions in the center-of mass frame, p(P A , S ⊥ ) + p(P B ) → J 1 (y c , P 1⊥ ) + J 2 (y d , P 2⊥ ) + X , where the polarized proton with the momentum P A and the transverse spin S ⊥ is moving in the +z-direction, while the unpolarized proton with the momentum P B is moving in the −z-direction, and we have the center-of-mass energy s = (P A + P B ) 2 . The produced two jets J 1 and J 2 have rapidities y c,d and transverse momenta P 1⊥ and P 2⊥ , respectively. These jets will be reconstructed via a suitable jet algorithm [59] and in the rest of the paper, we consider both of them to be anti-k T jets with jet radii R. In order to access the transverse motion of the partons inside the protons, we concentrate in the back-to-back region where the transverse momentum imbalance q ⊥ is small. Here we define the average transverse momentum P ⊥ of the two jets and the transverse momentum imbalance q ⊥ as follows where one has q ⊥ P ⊥ in the back-to-back region. The production of such back-to-back dijets is illustrated in Fig. 1. In the transversely polarized proton-proton collisions, the transverse spin vector S ⊥ of the incoming proton and the transverse momentum imbalance q ⊥ of the two jets will be correlated, as advocated in [42]. This correlation is accounted for in the Sivers function, which leads to a sin(φ q − φ S )-azimuthal modulation in the cross section between φ q and φ S , the azimuthal angles of q ⊥ and S ⊥ , respectively. Below we summarize the factorized formalisms for dijet production in both unpolarized and polarized proton-proton collisions, and we provide more details for the relevant ingredients in the next section. x v P a H S P J a P Z p J g E N G h 5 A P O q L G S 3 + h l N 9 N e u e J W 3 T n I K v F y U o E c j V 7 5 q 9 u P W R q h N E x Q r T u e m 5 g g o 8 p w J n B a 6 q Y a E 8 r G < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 k F M M 6 B h s a l m + a / l U S P i K i h Y n g 2 r j u t 1 P Y 2 N z a 3 i n u l v b 2 D w 6 P y s c n L R 2 n i q H P Y h G r T k g 1 C i 7 R N 9 w I 7 C Q K a R Q K b I e T u 7 n f f k K l e S w f z T T B I K I j y Y e c U W M l v 9 n P G r N + u e J W 3 Q X I O v F y U o E c z X 7 5 q z e I W R q h N E x Q r b u e m 5 g g o 8 p w J n B W 6 q U a E 8 o m d I R d S y W N U A f Z 4 t g Z u b D K g A x j Z U s a s l B / T 2 Q 0 0 n o a h b Y z o m a s V 7 2 5 + J / X T c 3 w J s i 4 T F K D k i 0 X D V N B T E z m n 5 M B V 8 i M m F p C m e L 2 V s L G V F F m b D 4 l G 4 K 3 + v I 6 a V 1 V v V r 1 9 q F W q T f y O I p w B u d w C R 5 c Q x 3 u o Q k + M O D w D K / w 5 k j n x X l 3 P p a t B S e f O Y U / c D 5 / A L C 7 j q E = < / l a t e x i t > P 1? < l a t e x i t s h a 1 _ b a s e 6 4 = " V 4 i E q j 5 u V r z G 4 L A T M r m 3 i / 4 D j N Q = " > A A A B + X i c b V B N S 8 N A E N 3 U r 1 q / o h 6 9 L B b B U 0 m k o N 6 K X j x W s B / Q h L D Z T t q l m 0 3 Y 3 R R K y D / x 4 k E R r / 4 T b / 4 b t 2 0 O 2 v p g 4 P H e D D P z w p Q z p R 3 n 2 6 p s b G 5 t 7 1 R 3 a 3 v 7 B 4 d H 9 v F J V y W Z p N C h C U 9 k P y Q K O B P Q 0 U x z 6 K c S S B x y 6 I W T + 7 n f m 4 J U L B F P e p a C H 5 O R Y B G j R B s p s G 1 v C j R v F 0 H u e i n I t A j s u t N w F s D r s e y t W K V M 6 f o D 6 z P H 9 x Y k 9 Q = < / l a t e x i t >P

2?
< l a t e x i t s h a 1 _ b a s e 6 4 = " P U C v v 5 L l B M c t l p U c i 2 M A a 9 C U y I s = " a t w 8 t O r t 2 z K O K j p D 5 + g S u e g K t d E 9 6 q A u o m i K n t E r e r N y 6 8 V 6 t z 6 W r R W r n D l F f 2 B 9 / g D d 4 p P V < / l a t e x i t > S ? < l a t e x i t s h a 1 _ b a s e 6 4 = " J c K j g h C l A 7 H a c b 6 u 0 t r 6 x u V X e r u z s 7 u 1 X 7 Y P D t k o y S a F F E 5 7 I b k g U c C a g p Z n m 0 E 0 l k D j k 0 A l H t z O / M w a p W C I e 9 S Q F P y Y D w S J G i T Z S Y F e 9 M d D 8 Y R r k X g o y n Q Z 2 z a k 7 c + B V 4 h a k h g o 0 A / v L 6 y c 0 i 0 F o y o l S P d d J t Z 8 T q R n l M K 1 4 m Y K U 0 B E Z Q M 9 Q Q W J Q f j 4 / f I p P j d L H U S J N C Y 3 n 6 u + J n M R K T e L Q d M Z E D 9 W y N x P / 8 3 q Z j q 7 8 n I k 0 0 y D o Y l G U c a w T P E s B 9 5 k E q v n E E E I l M 7 d i O i S S U G 2 y q p g Q 3 O W X V 0 n 7 v O 5 e 1 K / v L 2 q N m y K O M j p G J + g M u e g S N d A d a q I W o i h D z + g V v V l P 1 o v 1 b n 0 s W k t W M X O E / s D 6 / A F s u p O c < / l a t e x i t > J 1 < l a t e x i t s h a 1 _ b a s e 6 4 = " c X S A r 0 Y 2 G k t C 0 O l l 2 N d c u m 2 S L z k = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E U G 9 F L + K p o v 2 A N p T N d t I u 3 W z C 7 k Y o o T / B i w d F v P q L v P l v 3 L Y 5 a O u D g c d 7 M 8 z M C x L B t X H d b 6 e w s r q 2 v l H c L G 1 t 7 + z u l f c P m j p O F c M G i 0 W s 2 g H V K L j E h u F G Y D t R S K N A Y C s Y 3 U z 9 1 h M q z W P 5 a M Y J + h E d S B 5 y R o 2 V H u 5 6 X q 9 c c a v u D G S Z e D m p Q I 5 6 r / z V 7 c c s j V A a J q j W H c 9 N j J 9 R Z T g T O C l 1 U 4 0 J Z S M 6 w I 6 l k k a o / W x 2 6 o S c W K V P w l j Z k o b M 1 N 8 T G Y 2 0 H k e B 7 Y y o G e p F b y r + 5 3 V S E 1 7 6 G Z d J a l C y + a I w F c T E Z P o 3 6 X O F z I i x J Z Q p b m 8 l b E g V Z c a m U 7 I h e I s v L 5 P m W d U 7 r 1 7 d n 1 d q 1 3 k c R T i C Y z g F D y 6 g B r d Q h w Y w G M A z v M K b I 5 w X 5 9 3 5 m L c W n H z m E P 7 A + f w B y W 2 N f g = = < / l a t e x i t > J 2 < l a t e x i t s h a 1 _ b a s e 6 4 = " u R U 1 A t K c 7 X p M x T H N M y C s g 1 O r P 7 4 = " s C C D U T V K k + s W P t p l R q z g Q s C k 6 i I K Z s Q k f Q N z S k A S g 3 X R 6 9 w B d G G W I / k q Z C j Z f q 9 4 m U B k r N A 8 9 0 B l S P 1 W 8 v E / / y + o n 2 6 2 7 K w z j R E L L V I j 8 R W E c 4 S w A P u Q S m x d w Q y i Q 3 t 2 I 2 p p I y b X I q m B C + P s X / k 0 6 l T K r l x l 2 1 1 L x e x 5 F H Z + g c X S K C a q i J b l E L t R F D U / S A n t C z N b M e r R f r d d W a s 9 Y z p + g H r L d P o t C S s w = = < / l a t e x i t > FIG. 1. Illustration of back-to-back dijet production in transversely polarized proton-proton collisions: p(PA, S ⊥ ) + p(PB) → J1(yc, P 1⊥ ) + J2(y d , P 2⊥ ) + X. The polarized proton with momentum PA and transverse spin S ⊥ is moving in +z-direction, while the unpolarized proton with momentum PB is moving in −z-direction. We have jet rapidities y c,d and transverse momenta P 1⊥ and P 2⊥ , respectively. The dijet transverse momentum imbalance is defined as q ⊥ = P 1⊥ + P 2⊥ . Sivers asymmetry is generated due to the correlation between S ⊥ and q ⊥ .

A. Dijet unpolarized cross section
In the back-to-back region where q ⊥ P ⊥ , within the framework of soft-collinear effective theory (SCET) [60][61][62][63][64], one can write down a factorized form for the unpolarized differential cross section whereŝ = x a x b s is the partonic center-of-mass energy, N init is the corresponding spin-and color-averaged factor for each channel, while 1/(1 + δ cd ) arises from the symmetry factor due to identical partons in the final state. We have used the following short-hand notation In Eq.
(3), f unsub a (x a , k a⊥ , µ, ν) and f unsub b (x b , k b⊥ , µ, ν) are the so-called unsubtracted TMD PDFs, which carry the longitudinal momentum fractions x a,b and the transverse momenta k a⊥ and k b⊥ with respect to their corresponding proton. In our process, we have where y c , y d are the rapidities of the two leading jets. After performing Fourier transform for Eq. (3), we obtain the factorized formula in the coordinate b-space as follows dσ dy c dy d dP 2 Note that both the hard function H ab→cd and the global soft function S ab→cd are expressed in the matrix form in the color space and the trace Tr[· · · ] is over the color. Such factorization of the hard and soft function into matrix form is essential to capture evolution effects between the hard scale ∼ P ⊥ and the imbalance scale ∼ q ⊥ [65]. Here µ and ν denotes renormalization and rapidity scales, separately. The rapidity scale ν arises because both the TMD PDFs and the global soft functions have rapidity divergence [66,67], which are canceled between them as demonstrated below. This cancellation allows us to define rapidity divergence independentS ab→cd (b, µ) by where S ab (b, µ, ν) is the standard soft function appearing in usual Drell-Yan and SIDIS processes. This explicit redefinition allows us to subtract the rapidity divergence from the unsubtracted TMD PDFs to define the standard TMD PDFs f i (x i , b, µ) that are free of rapidity divergence as [68] f unsub Note that the properly-defined TMD PDFs f a (x a , b, µ) and f b (x b , b, µ) are no longer subject to the rapidity divergence and this is why there are no explicit ν-dependence in the arguments any more. Such properly-defined unpolarized TMD PDFs are the same as those probed in the standard SIDIS and Drell-Yan processes. The jet functions J c (P ⊥ R, µ) and J d (P ⊥ R, µ) in Eq. (6) describe the creation of anti-k T jets from the partons c and d, respectively. Finally, S cs c (k c⊥ , R, µ) and S cs d (k d⊥ , R, µ) are the collinear-soft functions. They describe soft gluon radiation with separations of order R along the jet direction, which can resolve the substructure of the jet. If one performs the integration over the azimuthal angle of the vector b, we obtain the following expression where J 0 is the Bessel function of order zero.

B. Dijet Sivers asymmetry
In the transversely polarized proton-proton collisions, the Sivers function will lead to a spin asymmetry in the cross section when one flips the transverse spin of the incoming proton. We thus define the difference in the cross section as d∆σ(S ⊥ ) = [dσ(S ⊥ ) − dσ(−S ⊥ )] /2. One can write down a similar factorized formula for such a spin-dependent differential cross section following Eq. (3), and it is given by where αβ is a two-dimensional asymmetric tensor with 12 = +1, and we have replaced the unpolarized TMD PDF in Eq. (3) by the Sivers function in the above equation following the so-called Trento convention [69], Note that we have also assumed that the global soft function S ab→cd (λ ⊥ , µ, ν) stays the same as that of the unpolarized collisions in Eq. (3). Although this is a reasonable assumption since the soft gluon radiation should be spin-independent [70,71], this has to be carefully checked. In fact, Ref. [72] shows in explicit calculations at one-loop level that soft functions in the polarized case can be different from the unpolarized counterpart beyond leading logarithmic accuracy, which is an indication of TMD factorization breaking. In this respect, our starting point Eq. (10) will be the best assumption at hand that takes a factorized form. We show the RG consistency for this factorized form, and we also demonstrate how we derive the process-dependent hard functions H Sivers ab→cd (P ⊥ , µ) for the polarized scattering. We leave a detailed study on the numerical impact of any TMD factorization breaking effects for future investigation.
Performing Fourier transform from the transverse momentum space into the b-space, we obtain where we have already used Eq. (7) to rewrite the unsubtracted unpolarized TMD PDF and Sivers function in terms of the properly defined versions which are free of rapidity divergence. Here f where we have used the fact that the integration in the first line would be proportional to b β , and we thus factored b β out explicitly in the second line 1 . The remaining part of the Sivers function is now denoted asf ⊥ a 1T (x a , b, µ). Note that for the same reason as explained below Eq. (8), we do not have the rapidity ν-dependence in the above equation. It is also instructive to emphasize thatf ⊥ a 1T (x a , b, µ) follows the same TMD evolution equations as the unpolarized TMD PDF f a (x a , b, µ), which enables us to evolve the Sivers function from some initial scale µ 0 to the relevant scale µ. On the other hand, at the initial scale µ 0 , the unpolarized TMD PDF f a (x a , b, µ 0 ) can be expanded in terms of the collinear PDFs f a (x a , µ 0 ). At a specific scale where the coefficient C a←i can be found in e.g. Refs. [68,73]. Likewise, Sivers functionf ⊥ a 1T (x a , b, µ) can be further matched onto the collinear twist-three Qiu-Sterman function T a,F (x 1 , x 2 , µ). At the scale µ b , one has the following expression for quark Sivers functionŝ where the matching coefficients at the NLO are given by [74][75][76][77][78] We now plug Eq. (13) into Eq. (12), and integrate over the azimuthal angle of the vector b, we obtain where J 1 is the Bessel function of order one, and we have used the identity withq ⊥ the unit vector along the direction of the imbalance q ⊥ . In general, the so-called single spin asymmetry (the Sivers asymmetry) A N for dijet production will be then given by Finally, since the Sivers function is not universal, one has to carefully include those non-universality or processdependence into the above formalism [43-45, 68, 79-81]. We have chosen to include all such process-dependence into the hard function H Sivers ab→cd (P ⊥ , µ), and this way the Sivers functions in Eq. (17) are the same as those probed in the SIDIS process. We explain in details how we derive the hard functions H Sivers ab→cd for different partonic processes in the next section.

C. Remarks
We will provide detailed expressions and discuss the evolution of all the relevant functions in the next section. Here, let us emphasize the following points on our factorized formalism: • Eqs. (6) and (12) are our proposed factorized formulas for dijet production in unpolarized and transversely polarized proton-proton collisions, respectively. They are the essential theoretical formalism we are using in the phenomenology section to compute the dijet Sivers asymmetry, which can be compared with the experimental data at the RHIC.
• It is important to emphasize that we have derived both Eqs. (6) and (12) within the SCET framework, in which the Glauber mode is absent. However, it is well-known that the inclusion of the Glauber modes will lead to factorization breaking. The factorization violation effects from Glauber gluon exchanging diagrams between two incoming nucleons have been discussed in [48,49,82,83]. In principle, such effects can be systematically accounted for in SCET by including explicitly the Glauber mode [84]. How exactly this works for dijet production remains to be investigated. In any case, the formalism we presented here would be a good starting point. This formalism incorporates the process dependence of the Sivers functions as outlined in [43-45, 79, 81], and also properly takes care of the QCD resummation and evolution effects. Thus in this formalism, we are able to study the energy and scale dependence of the Sivers asymmetry as measured in the experiment.
• There will be non-global structures from quantum correlations between in-jet and out-of-jet radiations: exclusive jet production will be sensitive on the correlation effects between in-jet and out-of-jet radiations, which is first discovered in [85]. The corresponding factorization and resummation formula involves multi-Wilson-line structures [86,87], which will give the non-linear evolution equation [88] for non-global logarithms (NGLs) resummation. The TMD factorization formula including such effects have been given in [56,89,90]. Numerically, the leading-logarithmic NGLs resummation can be solved using parton shower methods [85,[91][92][93] or BMS equations [94,95]. In our phenomenology, we have included the contributions from leading-logarithmic NGLs as discussed in Section V.
• Our formalism for unpolarized dijet production in Eqs. (6) is similar to those in [53,54]. Here, by taking the small-R limit, we refactorize the TMD R-dependence soft function [53,54] as the product of the R-independent global TMD soft function and the R-dependent collinear-soft function [55,56]. In addition, the R-dependent hard function in [53,54] has been further factorized into a R-independent hard function as above and the jet functions which naturally capture all the R-dependence. In this regard, the factorized formula presented here is more transparent and intuitive. Such refactorizations are essential to resum logarithms of R for small radius jets.
• After performing the refactorization mentioned in the above item, both the single logarithmic anomalous dimensions of the global and collinear-soft function not only depend on the magnitude | b| but also the azimuthal angle φ b of the vector b [55,56]. Especially, after taking into account QCD evolution effects the φ b integral is divergent in some phase space region. In order to regularize such divergences, we can first take φ b averaging in both the global and collinear-soft function, and then explicit φ b dependence will vanish. Therefore, one can avoid such divergence in the resummation formula directly. This φ b averaging method will not change the RG consistency at the one-loop order. The other methods to avoid such divergence have been discussed in [56], and no significant numerical differences are found at the NLL accuracy. The similar φ b averaging methods have also been used in [96][97][98] to simplify the calculation of the TMD soft function.

III. HARD FUNCTIONS IN UNPOLARIZED AND POLARIZED SCATTERING
In this section, we derive the hard functions for both unpolarized and polarized scatterings, i.e. H ab→cd (P ⊥ , µ) and H Sivers ab→cd (P ⊥ , µ) in Eqs. (9) and (17), respectively. They are matrices in the color space. We first review the results for the hard functions H ab→cd in the unpolarized scattering, which are well-known in the literature, see e.g. Refs. [70,99]. We then derive the hard function matrices H Sivers ab→cd in the polarized scattering case. These hard functions properly take into account the process-dependence of the Sivers functions [43-45, 68, 79-81]. To get started, we define the Mandelstam variables for the partonic scattering process, a(p 1 ) + b(p 2 ) → c(p 3 ) + d(p 4 ), as followŝ Unpolarized scattering amplitudes for the qq → qq subprocess. From left the right, the scattering amplitude is provided for the t-and u-channel processes.
where ∆y = y c − y d is the rapidity difference of the two jets. In the following, the expressions for the hard functions will be written in terms of these Mandelstam variables.
A. Unpolarized Hard Matrices

Four quark subprocesses
We start with the partonic subprocesses that involve four quarks, such as qq → qq. In Tab. I, we organize each of the four quark subprocesses into a color basis. The color basis operators acting on particles i and j are denoted as Γ n,ij which are used to generate the hard and soft matrices. For the four quark interactions, two operators, n = 1, 2, are required to span the color space. As seen in the table, this results in 12 total color matrices. Using the fact that hard function for the unpolarized case is invariant under the charge conjugation, the bottom row can easily be computed from the top row. Furthermore, once the hard matrices have been calculated for the first column, crossing symmetry can be applied in order to obtain the hard color matrices for the second and third column. It is then only necessary to explicitly calculate the hard matrices for the subprocesses associated with the color basis Γ n,31 Γ n,42 . For our calculation, we follow the conventions used in Refs. [70,99] to choose Γ 1,ij = (t a ) ij and Γ 2,ij = δ ij , so that the color basis is spanned by the orthogonal basis We note that other bases have been used in the literature [100]. We now explicitly perform the calculation for the qq → qq , qq → q q, and qq → qq subprocesses. For these subprocesses, we can write where we have suppressed the ab → cd subprocess label. The subscript in the M terms denotes the relevant Mandelstam variable (t orû) for the channel that contributes to the subprocess as shown in the Fig. 2. To arrive at this expressions, we have separated the color parts from the kinematic parts (denoted with the superscript kin). These kinematic scattering amplitudes are defined by We can now decompose these scattering amplitudes in color space as where To obtain the expressions in Eq. (27), we have exploited the orthogonality of our chosen color basis in Eqs. (21) and (22). Then we will have |M| 2 as where the hard matrix is given by and the leading order soft matrix as The hard matrices of the four quark processes in Γ 31 Γ 42 color basis in Tab. I are given by We find these results to be consistent with the expressions in [70]. The remaining hard functions can be obtained from crossing symmetries.

Two quarks and two gluon subprocesses
In Tab. II, we provide a list of subprocesses involving two quarks and two gluons with the color basis operators Γ ab n,ij . For the two quark and two gluon interactions, three operators, n = 1, 2, 3, are required to span the color space. A convenient choice for the computation is the set of orthogonal operators (primed), which has the corresponding orthogonal basis, At the same time, we find that the final expressions for the hard matrices take a simpler form when one uses the non-orthogonal basis used in Refs. [70,99,100] by defining the basis operators to be (unprimed) The corresponding basis is given by We note that the normalization of θ 3 in [100] differs from the normalization of Refs. [70,99] by a factor of 2. For the choice of basis in Eq. (37), the LO soft matrix is given by In order to exploit the orthogonality condition of the primed basis in Eq. (35), but still provide a simple expression for the hard matrices using the unprimed basis in Eq. (37), we first compute the hard matrices in the primed basis then obtain the results in the unprimed basis using the relation We now perform the calculation for the hard matrices for the qq → gg process in the primed orthogonal basis. The scattering amplitude for this subprocess can be written in color space as where The hard matrix in the primed basis can therefore be computed as Finally, we now use Eq. (39) to obtain the simplified hard functions in the unprimed basis as The hard matrices for other subprocesses involving two quarks and two gluons, such as qg → qg, can be obtained from this expression using crossing symmetries.

Four gluon subprocesses
For the four gluon subprocesses, gg → gg, we follow the work in Refs. [70,99] to use the following over-complete basis We note that a six dimensional basis was chosen in [100]. Using this basis in Eq. (45), one can show that the hard matrix takes the following form The LO soft matrix for this channel is given in Appendix C of [99] for this basis as

B. Polarized Hard Matrices
As we have emphasized in the previous section, Sivers function is non-universal. The well-known example is the sign change between the Sivers function probed in SIDIS and that in Drell-Yan (DY) process [32][33][34], Such a sign change can be easily taken care of in describing the Drell-Yan Sivers asymmetry, where H(Q, µ) is the hard function in the Drell-Yan process, and we have applied Eq. (48) in the second step. In other words, if we use the SIDIS Sivers function in a Drell-Yan process, we shift the minus sign (or the process-dependence) into the hard function. For the partonic subprocesses in the hadronic dijet production, one has much more complicated process-dependence for the Sivers functions involved. This can be seen from the highly nontrivial gauge link structure which has been derived in [80] in the definition of the TMD PDFs. Even in these complicated processes, one can incorporate such process-dependence of the Sivers functions into modified hard functions as in Eq. (49) [43-45, 79, 81]. We follow a similar procedure in this section to include this process-dependence of the Sivers functions into the hard functions in the matrix form.
In Fig. 3, we demonstrate the factorization between the Sivers function and modified hard functions. Unlike the unpolarized case, the contributions of the Sivers asymmetry are given by considering the attachment of an additional collinear (to the incoming hadron) gluon to three of the external legs. Such a gluon is part of the gauge link in the definition of the Sivers function, and it is the imaginary part of the Feynman diagram (related to the so-called soft gluonic pole) that contributes to the process-dependence of the Sivers function.
It is important to note that the additional gluon leads to additional complications so that naive crossing symmetry cannot be used to relate one hard function to another, as in the unpolarized case studied above. These complications occur because the contributions to the Sivers asymmetry are only given by attaching the additional gluon to three of the four external legs. Furthermore, since the sign of the interaction (imaginary part) with the external gluon is opposite for quarks and anti-quarks, this sign must also be accounted for when applying crossing symmetry or charge conjugation.

Four quark subprocesses
As in the unpolarized case, the bases for four quark subprocesses are given in Tab. I. As discussed above, one cannot naively apply crossing symmetry to obtain hard matrices of a general polarized subprocess. For the polarized four quark subprocesses, however, only the sign of each color factor changes under charge conjugation. Therefore, the hard matrices for the bottom row of Tab. I can be obtained from the results from the top row of this table with the addition of a minus sign. To demonstrate how H Sivers ab→cd are derived, we explicitly perform the calculation for the qq → qq , qq → q q, and qq → qq subprocesses as we did for the unpolarized case. Afterwards, we provide the expressions for the remaining subprocesses. To start, it is important to remind ourselves that a non-vanishing Sivers asymmetry requires initial/final state interactions generating a phase. Because all initial and final partonic states relevant for dijet production are colored, both initial and final state interactions have to be taken into account. Such interactions would generate non-trivial gauge link structures, see e.g. Refs. [44,80,101]. On the left side of Fig. 3, as an example, we show all possible diagrams with one gluon exchange between the remnant of the polarized proton and the qq → qq hard scattering part, which contribute to the Sivers asymmetry. Now with the presence of the extra gluon scattering (first order of the gauge link expansion), the diagram at the left side of the cut will be denoted as M Sivers,a j , while the right side is same as the unpolarized case denoted as M † . Here a is the color for the attached gluon, j is the color index for the incoming quark with momentum x a P A on the left side of cut line, while the color index for the incoming quark on the right side of the cut line is given by 1 like in the previous section. In contrast to the unpolarized correlation function, quarks j and 1 do not need to have the same color, because of the presence of the gluon from the gauge link. Now we perform the following expansion to obtain the hard matrix |M Sivers | 2 for the polarized case, where t a 1j will be included into the quark-quark correlator in the polarized proton to become ∼ P S|ψ 1 n·A a t a 1j ψ j |P S , see e.g. Ref. [80,81,102]. From Eq. (50), we thus derive At the same time, we use the convention that N init in the polarized and unpolarized cases are the same. Therefore, the factor of 1/N c in Eq. (51) is absorbed into N init . With that in mind, to arrive at the correct normalization of the polarized hard function, we thus obtain which is demonstrated on the right-hand side of Fig. 3. Now we need to project M Sivers,a j and M † into the color basis separately. The polarized scattering amplitude M Sivers,a j can be written as where M kin t and M kin u are the same as the expressions in Eqs. (24) and (25). From left to right on the top line of this expression, these terms give the scattering amplitudes for the initial-state, final-state 1, and final-state 2 interaction for the t-channel, corresponding to the first three diagrams of Fig. 4 in the same order. Likewise from left to right on the bottom line, the terms give the scattering amplitude for the initial-state, final-state 1, and final-state 2 interaction for the u-channel, corresponding to the last three diagrams of Fig. 4 in the same order. Using the Feynman rules for the gauge link color factors given in Fig. 6 of [81], we easily arrive at Eq. (53) from these diagrams. From the unpolarized scattering amplitude given in Eq. (23), we write the conjugate amplitude as Analogous to the unpolarized scattering amplitude, the scattering amplitude can be decomposed into the orthogonal basis given in Eq. (21) as where we have After performing this decomposition, we can now write where H Sivers ab→cd is given by and S is the same as the unpolarized case. From these expressions, we can obtain the polarized hard matrices for the qq → qq , qq → q q, and qq → qq subprocesses as Since qq → qq subprocess receives contributions from both t-and u-channels (as well as their interference), its expression is the most complicated among the three subprocesses computed. One can show that after performing the trace with the soft color matrix, the expressions are consistent with the squared amplitude of [81]. The color matrices for the remaining four quark subprocesses in the top row of Tab. I can be computed in the same spirit and we obtain the following expressions After performing charge conjugation, the hard color matrices for the subprocesses in the bottom row of Tab. I can be obtained from these expressions.

Two quarks and two gluon subprocesses
All twelve of the two quark and two gluon subprocesses are given in Tab. II. As we have mentioned in Sec. I, we neglect the gluon Sivers contribution in this paper. This means that all subprocesses with a gluon incoming from the polarized proton will be neglected. There are then six remaining subprocesses to compute. However, we find that under charge conjugation, the polarized hard functions once again only change by an overall minus sign. Thus, we only need to perform the calculation for three of the hard matrices.
In order to further demonstrate our method for calculating the polarized hard matrices, we now perform the calculation for the qq → gg subprocess. We then provide the expressions for the remaining hard matrices. For the unpolarized process the scattering amplitude has three channels. After the addition of the external gluon, there are then nine polarized process to be considered. At the cross section level, this results in 27 hard interactions which need to be considered. Despite this complication, we can once again write where M kin = M kin s + M kin t + M kin u in general receives contribution from different channels as above. The kinematic parts can be trivially extracted by Therefore the unpolarized hard matrices can be constructed simply by 2 In these expressions, we have suppressed the subprocess subscript since these expressions are true for all subprocesses with a one-dimensional color space. The differential cross section is then given by where in the second line we have defined C u = Tr θ 1 θ † 1 . Similarly, for the polarized hard matrix, we can write where C Sivers M kin = Tr M Sivers,a t a j1 θ † 1 . Therefore, the hard functions of the polarized and unpolarized scatterings are related by an overall color constant, Here, C Sivers can further be decomposed into color factors arising from gauge link gluons interacting with different external colored partons, as seen in [81,[103][104][105].

C. Evolution equations
Hard functions can be related to the Wilson coefficients C Γ I in the color basis {θ I } of section III by H IJ = Γ C Γ I C Γ * J . Here Γ represents different helicity states of the incoming and outgoing particles. Explicit expressions of the Wilson coefficients at next-to-leading order can be found in [70,99], but we do not present them as we are only using the tree-level hard functions for our study. We do, however, include the renormalization group (RG) evolution of the hard functions coming from the 1-loop anomalous dimensions. Then the Wilson coefficients satisfy the RG evolution equations [70,99,106,107] µ d dµ Here, γ cusp = αs π + · · · is the cusp anomalous dimensions and c H = C a + C b + C c + C d . The non-cusp anomalous dimension is defined as where γ i µ [α s (µ)] = αs π γ i + · · · , with γ q = 3 2 C F and γ g = β0 2 . Lastly, the matrix M takes the form where s 12 = s 34 =ŝ, s 13 = s 24 =t, and s 14 = s 23 =û and From the RG evolution of the Wilson coefficients given in Eq. (87), we can arrive at the RG evolution equations for hard matrix H as where Γ H is given by

IV. QCD RESUMMATION AND EVOLUTION FORMALISM
In this section, we present the renormalization group (RG) equations for the rest of the key ingredients in the factorized formalism. These include the TMD PDFs, global soft functions, jet functions, and collinear-soft functions. After presenting their NLO perturbative results and RG evolution equations, we check the RG consistency. In the end, we present our resummation formula for dijet production.

A. TMDs and global soft functions
The unsubtracted TMD PDFs in the factorized formula in Eq. (6) describe the radiation along the incoming beams. They satisfy the RG evolution equations where its µ-and ν-anomalous dimensions are given by As we will see in this subsection, the rapidity divergences of the unsubtracted TMDs will be exactly canceled by the rapidity divergences of the global soft functions, which will allow us to identify the standard TMDs with subtracted rapidity divergence as in Eq. (8) above. Suppressing the label ab → cd for convenience, the global soft functions up to 1-loop are given by where [108] I The explicit matrix forms of tree-level soft functions in Eq. (97) for some color basis {θ I } can be computed as which is equivalent to the matrix forms of the LO soft functions found in section III. The matrix T i · T j of the eq. (98) was also computed in the color bases used in section III and can be found in [70,99]. The renormalized global soft functions satisfy the RG evolution equations From Eqs. (97) -(102) and using i T i = 0, we then find where M was given in Eq. (89) and we promoted αs π → γ cusp , which is consistent with the factorization consistency relation below. Note that Eq. (107) is strictly real and the imaginary term ∼ iπ cancels exactly with the imaginary term found in M .
We note that Γ S ν ∼ I and that this is expected as the hard functions do not have any rapidity divergence. Thus, we can write which has the same rapidity anomalous dimensions as the back-to-back soft functions S ab (b, µ, ν) found in standard Drell-Yan and SIDIS process [67]. As expected, the rapidity divergence of the global soft function S(b, µ, ν) in Eq. (109) exactly cancels the rapidity anomalous dimensions for the unsubtracted TMDs f a (b, µ, ν) and f b (b, µ, ν) given in Eq. (96). Therefore, as discussed in the introduction, we can defineS(b, µ) absent of the rapidity divergence such that Then as in Eq. (8), S ab (b, µ, ν) is combined with the unsubtracted TMDs to identify standard TMDs free of the rapidity divergences.

B. Jet and collinear-soft functions
Both jet and collinear-soft functions describe the radiation which resolves the produced jets. The jet functions [109,110] encode the collinear radiations inside anti-k T jet with radius R. The NLO expressions are given by where the algorithmic dependent terms d i for anti-k T algorithm are The jet functions satisfy the RG evolution equations where the anomalous dimension is given by The collinear-soft functions [55,56] describe the soft radiation along the jet direction and resolves the jet cone R. The NLO expressions are given by The collinear-soft functions satisfy the RG evolution equations where its anomalous dimension takes the form C. RG consistency at 1-loop With the anomalous dimensions presented for all the ingredients, we now show that our factorized formula given in Eq. (6) satisfy the consistency relations for the RG evolutions. The cancellation of the rapidity divergences was already checked around Eq. (109). We also expect µ-divergence of the various functions to cancel and satisfy the consistency equation From Eqs. (91), (92), (104), (107), we immediately find at 1-loop, One can then easily check from the µ-anomalous dimensions of the other functions given in Eqs. (95), (115), (118) that Eq. (119) is explicitly satisfied at 1-loop.

D. Resummation formula
Based on the above discussions and RG renormalization group methods in SCET, we can now derive the expression for the all-order resummed result. Explicitly, we calculate the cross section at the NLL accuracy, where we will use the two-loop cusp and one-loop single logarithmic anomalous dimension and the matching coefficients are kept at leading order. On the other hand, the color structures inside the hard and soft function will mix with each other under the RG evolution, which was first studied in [65]. In this paper, we will apply the same methods in [70] to solve the RG equations. For the unpolarized cross section, the resummation formula has the form as follows: where λ K is the eigenvalue of the matrix M IJ in the hard anomalous dimension (87) and H KK andS K K are the hard and soft function in the diagonal basis as defined in [70]. In our numerical calculation, we use the LAPACK library [111] to obtain their value at different phase-space points. We have applied the b * -prescription to prevent the Landau pole from being reached in the b-integral. Here, we define b * as where b max is chosen [112] to be 1.5 GeV −1 . The nonperturbative Sudakov factor in Eq. (121) was fitted to experimental data in [113]. The extracted functions are given by We also incorporate NGLs resummation effects included by the function U c,d NG . In order to include NGLs resummation effects at NLL accuracy, we also need to consider the extra one-loop single logarithmic anomalous dimensionΓ from the non-linear evolution parts. However, in [86,87] this anomalous dimension was shown to cancel between the jet and collinear-soft function up to two-loop order. The explicit operator-based derivation of RG consistency includinĝ Γ can be found in [56,90,114]. In the large N c limit, the non-linear evolution equation can be solved using the parton shower algorithm [115]. Especially, at the NLL accuracy the evolution is totally determined by the one-loop anomalous dimensionΓ, which is equivalent to the one appearing in the light jet mass distribution at the e + e − collider. Therefore, we can use the same fitting function form given in [85] to capture NGLs resummation contributions after setting proper initial and final evolution scales. In our case, these two scales are the jet scale µ j and the collinear-soft scale µ cs . Explicitly, the function is where the superscript k = q and g denote the (anti-)quark and gluon jet, respectively, and with C q = C F and C g = C A . The parameters a, b and c are fitting parameters which are given as a = 0.85 C A , b = 0.86 C A and c = 1.33. The variable u = 1 β0 log αs(µcs) αs(µj ) is the evolution scale measuring the separation of the scales µ cs and µ j . As we have done for the unpolarized cross section, we also derive a similar resummation formula for the spindependent cross section where at the NLL accuracy we keep the LO matching coefficient in Eq. (16). It involves the parametrization for the Sivers function, which depends on the collinear Qiu-Sterman function T q,F (x a , x a , µ b * ) and a different non-perturbative Sudakov factor S s NP . The relevant parametrization has been determined from a recent global analysis of the Sivers asymmetry of SIDIS and Drell-Yan processes [116]. The non-perturbative Sudakov factor is given by

V. PHENOMENOLOGY
In this section we will present the numerical results using the resummation formula in Eqs. (121) and (125), where intrinsic scales for the hard, jet and collinear-soft function are chosen as In the numerical study, we will focus on the Sivers asymmetry for the dijet production at the RHIC with √ s = 200 GeV, where the jet events are reconstructed by using anti-k T algorithm with jet radius R = 0.6. The transverse momentum P ⊥ and the rapidity y c,d of jets are For the unpolarized proton, we use the HERAPDF20NLO parton distribution functions [117]. The numerical Bessel transforms in Eqs. (121) and (125) are performed using the algorithm in [118]. Furthermore, the Eq. (9) is derived after neglecting the power corrections from O(q 2 ⊥ /P 2 ⊥ ). In other words, in the large q ⊥ region, the full results should include corrections from the so-called Y -term, which can be obtained from perturbative QCD calculations [119]. In this paper we focus on the contribution from back-to-back dijet production. In order to select such kinematics, we require the transverse momentum q ⊥ for the dijet system |q ⊥ | < q cut ⊥ . In the numerical calculations, we fix the value of q cut ⊥ = 2 GeV. As shown in the Fig. 1, the transverse-polarized proton moves on +z-direction and its spin points to +y-direction with φ S = π/2. The transverse momentum vector q ⊥ lies in the x − y plane, and the Sivers asymmetry is defined as the difference of the events between q ⊥,x > 0 and q ⊥,x < 0 hemispheres, that is the same as the measurements by STAR collaboration [41]. Explicitly, we have with dPS = dy c dy d dP ⊥ δ(y sum − y c − y d ) represents the transverse momenta and rapidities integral for dijets. In the numerator, the φ q -integral with θ(cosφ q ) and θ(−cosφ q ) corresponds q ⊥,x > 0 and q ⊥,x < 0, respectively. In the Fig. 5, we show the numerical results of the Sivers asymmetry for dijet processes, where we neglect the charm and bottom jet events. The red and blue curves represent the asymmetry contributed from u-and d-quark Sivers function, respectively. As is expected, we find that the asymmetry is enhanced in the large y sum region, i.e. the forward scattering region, due to the larger fractional contribution of Sivers function in the valence region. Besides, the contributions from u-and d-quark Sivers function are opposite from each other, which causes a huge cancellation of the asymmetry, as shown by the black curves in Fig. 5.
In the calculation, most of the asymmetries come from the partonic scattering process qg → qg where the initial quark comes from the polarized proton. Especially, the more forward jet is associated with the parton from the polarized proton moving in the same direction. Hence, if we can tag parton species initiating the more forward jet, then we can separate u-and d-quark Sivers functions and avoid the accidental cancellation as shown in the left plot of Fig. 5.
In order to achieve jet flavor separation mentioned above, one possible method is applying the electric charge information of jets, which has been proposed in [50,58,120]. In this paper, we will use the standard jet electric charge definition given in [121,122] where z h is the transverse momentum ratio between hadrons and the jet. κ is an input parameter, which is fixed by κ = 0.3 [58] in our calculations. As shown in [58], after measuring the jet charge information, the theory formula is slightly modified by replacing the jet function J i (P ⊥ R, µ) in Eq. (17) by the charge-tagged jet function G i (Q κ , P ⊥ R, µ) as with the normalization as ∞ −∞ dQ κ G i (Q κ , P ⊥ R, µ) = J i (P ⊥ R, µ) required by the probability conservation. Here we only replace the more forward jet function with the charge-tagged jet function, which corresponds to the insertion of the step function. We define the jet charge bin fraction as Then the Sivers asymmetry A N in different jet charge bins is given as, in terms of jet charge bin fraction where we suppress the phase space integral shown in Eq. (129). The index i denotes the parton species initiating the more forward jet. Here we use the same jet charge bins defined in [58], where +, − and 0 indicate Q κ > 0.25, Q κ < −0.25 and |Q κ | < 0.25 bins, separately. Such jet charge bin fraction can be fitted from the unpolarized cross section for back-to-back dijet events at the RHIC. In [50], the authors have shown the preliminary results from the measurements as κ = 0. In the theory calculation, one can use Monte-Carlo event generators such as Pythia8 [123] to estimate these numbers. In the Tab. III we give the results of jet charge bin fractions r ±,0 i for various jet flavors used in our numerical calculations, where the jet charges are defined using all charged hadrons inside the jet. In the right plot of Fig. 5 we show the result of A N within the different jet charge bins. After selecting the charge of the more forward jet Q κ > 0.25, the contribution from the u-quark Sivers function is enhanced compared to the case without the jet charge measurement (the black curve in the left plot). A similar size enhancement from the d-quark Sivers function is also observed in Q κ < −0.25 charge bin as shown by the blue curve. Besides, we find the Sivers asymmetries from Q κ > 0.25 bins are positive and Q κ < −0.25 bins are negative, which are consistent with the preliminary STAR measurements [50]. In the forward region, the Sivers asymmetry can achieve O(0.01%), and size of our calculation is also around the same order of the data. Taken together, our calculation suggests that the dijet production at the hadron collider is an important process to extract the information about the Siver function and deserves further studies on the theoretical framework about the remarks discussed in II C.

VI. CONCLUSIONS
We study the single spin asymmetries of dijet production in the back-to-back region in transversely polarized protonproton collisions. In the back-to-back region, the dijet transverse momentum imbalance q ⊥ is much smaller than the transverse momentum P ⊥ of the jets. In this case, the conventional perturbative QCD calculations in the expansion of coupling constant α s generate large logarithms in the form of α n s ln m P 2 ⊥ /q 2 ⊥ with m ≤ 2n − 1, which have to be resummed in order to render the convergence of the perturbative computations. We propose a QCD formalism in terms of transverse momentum dependent (TMD) parton distribution functions for dijet production in both unpolarized and polarized proton-proton collisions. Such a formalism allows us to resum the aforementioned large logarithms, and further takes into account the non-universality or process-dependence of the Sivers functions in the case of the transversely polarized scattering. It is well-known that hadronic dijet production in back-to-back region suffers from TMD factorization breaking effects. Thus, to write down the QCD "seemingly factorized" formalism for resumming large logarithms mentioned above, we make a couple of approximations. First of all, we neglect the Glauber mode in the formalism which are known to be the main reason for the TMD factorization breaking. Secondly, we have assumed that the soft gluon radiation that is encoded in the global soft function in our formalism is spin-independent, i.e., they are the same between the unpolarized and polarized scatterings. Since the precise method for dealing with the TMD factorization breaking effects is still not known, we feel that the proposed formalism in this paper is a reasonable starting point for further investigation.
With such a formalism at hand, we compute the Sivers asymmetry for the dijet production in the kinematic region that is relevant to the proton-proton collisions at the Relativistic Heavy Ion Collider (RHIC), and find that the spin asymmetry is very small due to the cancellation between u-and d-quark Sivers functions, which are similar in size but opposite in sign. However, we find that the individual contribution from u-and d-quark Sivers functions can lead to an asymmetry of size O(±0.05%) in the forward rapidity region, which seems feasible at the RHIC. Motivated by this, we compute the Sivers asymmetry of dijet production in the positive and negative jet charge bins, i.e., when the jet charge Q κ for the jet with the larger rapidity of two is in the bins Q κ > 0.25 and Q κ < −0.25, respectively. By selecting the positive (negative) jet charge bin, we enhance the contribution from u-(d)-quark Sivers function and thus enhance the size of the asymmetry. Our calculation shows that Sivers asymmetries in such positive (negative) jet charge bins lead to asymmetries of size O(+0.01%) (O(−0.01%)), respectively. The sign of such asymmetries seem to be consistent with the preliminary STAR measurements at the RHIC. The size of our calculations is also around the same order of the experimental data. This give us a great hope to further investigate the single spin asymmetries for hadronic dijet production at the RHIC. Note added: While this work was being written up, we noticed a similar work [72] appears on arXiv. The authors investigate process dependent factorization violation from the soft gluon radiation. Their method is different from our approach. We assume a factorized form for the spin-dependent cross section, which we demonstrate to be renormalization group consistent. Within this factorized form, we explicitly calculate the process dependent polarized hard function in the matrix form. Besides, in the numerical calculations we include quark Sivers functions in all the partonic channels. We believe these two studies are complementary with each other.