Azimuthal Angular Asymmetry of Soft Gluon Radiation in Jet Production

We investigate the impact of soft gluon resummation on the azimuthal angle correlation between the total and relative momenta of two energetic final state particles (jets). We show that the initial and final state radiations induce sizable $\cos(\phi)$ and $\cos (2\phi)$ asymmetries in single jet and dijet events, respectively. We numerically evaluate the magnitude of these asymmetries for a number of processes in collider experiments, including diffractive dijet and dilepton production in ultraperipheral $pA$ and $AA$ collisions, inclusive and diffractive dijet production at the EIC and inclusive dijet production in $pp$ collisions at the LHC. In particular, the $\cos (2\phi)$ asymmetry of perturbative origin can dominate over the primordial asymmetry due to the linearly polarized gluon distribution.


I. INTRODUCTION
Jet processes are the most abundant events in hard hadronic collisions and have been under intensive investigations at various colliders, see, e.g., Refs. [1][2][3][4][5][6][7][8]. Among them, the productions of dijet and jet plus a color-neutral particle (such as the Higgs boson) are characterized by distinct final states where two energetic particles (jets) are almost back-to-back in the transverse plane perpendicular to the beam direction. Deviations from the exactly back-to-back configuration are in general expected. In other words, the total transverse momentum of the two outgoing systems q ⊥ = k 1⊥ + k 2⊥ is typically small but nonzero, as illustrated in Fig. 1 for the dijet case. The cross section then depends on the angle φ between q ⊥ and the dijet relative momentum P ⊥ = ( k 1⊥ − k 2⊥ )/2 dσ dP ⊥ dq ⊥ dφ = σ 0 + cos(φ)σ 1 + cos(2φ)σ 2 + · · · . (1) The coefficients σ 1 , σ 2 , · · · often encode novel partonic structures of the target that are important in the study of nucleon tomography at the future electron-ion collider (EIC) [9][10][11]. A primary example is exclusive diffractive dijet production in γ ( * ) p scattering where q ⊥ is provided by the recoil momentum of the target. It has been predicted that the 'elliptic' gluon Wigner distribution generates a cos(2φ) asymmetry [12][13][14][15][16][17]. Another example is the inclusive dijet production in DIS, where q ⊥ comes from the intrinsic transverse momentum of gluons in the target. In this case, the so-called linearly polarized gluon distribution can generate a cos(2φ) asymmetry in the dijet system [18][19][20][21][22][23][24]. However, the momentum imbalance q ⊥ can simply come from perturbative initial and final state radiations which have nothing to do with nontrivial parton distributions inside the target. Depending on kinematics, this can affect or even dominate the coefficients σ 1 , σ 2 , · · · when P ⊥ q ⊥ . The reason is that the radia- 1. Dijet in transverse plane perpendicular to the beam direction at hadron colliders. The dijet total transverse momentum q ⊥ = k 1⊥ + k 2⊥ , which is due to multiple soft gluon emissions, is much smaller than the individual jet momentum P ⊥ = ( k 1⊥ − k 2⊥ )/2. tive corrections are enhanced by large double logarithms (α s ln 2 P 2 ⊥ /q 2 ⊥ ) n . While the resummation of these logarithms is well understood for the angular independent part σ 0 [25][26][27][28][29][30][31][32][33][34][35], that for the angular dependent part has been discussed much less frequently in the literature. In a series of papers by Catani et al. [36,37], it has been observed that the resummation for σ 1 , σ 2 , · · · can be done in the Fourier space q ⊥ → b ⊥ using the same Sudakov factor as for σ 0 . An interesting new feature is that although the angular dependent cross section is singular 1/q 2 ⊥ in fixed-order calculations with no compensating virtual correction, the resummed cross sections σ 1,2 are well-behaved as q ⊥ → 0. The goal of this paper is to study this resummation in detail and make quantitative predictions for azimuthal asymmetries cos(nφ) that can be compared with the existing and future experimental data. We shall consider a variety of processes, including dijet productions in diffractive and inclusive processes, and jet plus color-neutral particle production. A brief summary of our results has been published in Ref. [38].
An important feature of the final state radiation is that the emitted soft gluons tend to be aligned with jet direc- Since k i⊥ 's tend to point to jet directions, the same is true for q ⊥ , resulting in a sizable anisotropy cos(2φ) .
tions, see Fig. 2. Those emitted inside jet cones become part of the jets, so one needs to carefully treat emissions slightly outside the jet cones. Since q ⊥ is the recoil momentum against these gluons, it also points towards jet directions on average. This naturally generates a positive cos(2φ) asymmetry of purely perturbative origin in dijet events [38]. A recent measurement by the CMS collaboration [39] indicates that the magnitude of this effect is sizable. Depending on kinematics, it can completely overshadow the 'intrinsic' azimuthal correlations generated by nontrivial parton distributions. For the search of the elliptic Wigner distribution in exclusive dijet production mentioned above, one can avoid this problem by measuring the correlation between P ⊥ and the nucleon recoil momentum (instead of q ⊥ ) as originally suggested in [12]. However, in inclusive dijet production, it is not possible to cleanly separate the contribution from the linearly polarized gluon distribution, unless one has an accurate control of the perturbative backgrounds.
Our discussions in this paper are connected to other recent developments in the field. The correlation between q ⊥ and P ⊥ also measures the correlation between final state jets. The combination of this study with three particle correlations in the final state recently proposed in Refs. [40][41][42] shall open a new avenue to study the QCD dynamics of gluon radiation. In this regard, the nonglobal logarithms (NGLs) [26,43] can also contribute to the observables we consider, although their numerical impact might be limited for the relevant kinematics [31]. More broadly, the perturbative contribution to the cos(2φ) azimuthal asymmetries has been studied for various processes [36][37][38][44][45][46][47][48][49]. In particular, it may shed light on the QCD factorization and resummation for power corrections in hard scattering processes [47,[50][51][52][53].
This paper is organized as follows. In Section II, we consider processes in which the final state consists of one jet and one color-neutral particle. In subsection IIA, we consider lepton-jet correlations in ep scattering in the laboratory frame. In subsection IIB, we consider photonjet production in pp collisions. The former can be studied at the EIC, whereas the latter can be studied at RHIC and the LHC. Since there is only one jet in the final state, we expect that the dominant asymmetry is of the form cos(φ).
In Section III, we study dijet production. In subsections IIIA and IIIB, we consider diffractive and inclusive dijet photo-production processes, respectively. As mentioned above, the dominant asymmetry is cos(2φ) in this case. Then in subsection IIIC, we consider inclusive dijet production in pp collisions specifically focusing on the (most complicated) gg → gg channel.
In Section IV we give a detailed analysis of dijet electro-production in DIS. When the photon is virtual, the linearly polarized gluon distribution gives an additional contribution to the cos(2φ) asymmetry. We numerically compare the respective contributions to the asymmetry from nonperturbative and perturbative mechanisms.

II. JET PLUS COLOR-NEUTRAL PARTICLE IN THE FINAL STATE
In this section, we discuss the final states with a jet and a color-neutral particle. The soft gluon radiation comes only from the jet and the incoming parton(s). The dominant azimuthal asymmetry is then the cos(φ) term. We will first study lepton plus jet production in ep collisions and then extend to photon plus jet production in pp collisions. Similar studies can also be carried out for jet plus Higgs boson or Z/W boson productions at the LHC.

A. Lepton and Jet Correlation in ep Collisions
At leading order a lepton scatters off a quark through virtual photon exchange in t-channel and produces a quark jet in the final state e(k) + q(p 1 ) → e (k ) + jet(k J ) + X . ( In the laboratory frame, the final state lepton and jet are mainly back-to-back in the transverse plane perpendicular to the beam direction. This process has recently attracted significant interest because it can provide a novel way to study the transverse momentum dependent (TMD) quark distribution in the nucleon [32,33,[85][86][87][88][89]. It has also motivated experimental efforts to re-analyze the existing HERA data [90,91]. The virtuality of the photon defines the hard-scattering process. To leading order, the differential cross section can be written as where , y is the rapidity of the final state lepton in the laboratory frame. Following the notations in Introduction, we have define the difference and total transverse momenta for the two final state particles: [Below we often omit an arrow on transverse vectors.] In the above equation, x represents the momentum fraction of the incoming nucleon carried by the quark, f q (x) for the quark distribution function. The Mandelstam variablesŝ,t and u are defined as usual for the partonic sub-process, in At one-loop order, q ⊥ can be nonzero due to the emission of a soft gluon with momentum k ⊥g [33]. Integration over the phase space of the emitted gluon is explained in Appendix A. The result is where and φ is the azimuthal angle between q ⊥ and P ⊥ . Note that the coefficients c n in general depend on q ⊥ . But the dependence is power-suppressed where the integer c (usually c = 1 or 2) depends on both n and the process under consideration. In (4), we recognize at least two sources of such power corrections. First, when q ⊥ is small but nonvanishing, the soft gluon rapidity is subject to kinematical constraints y min < y g < y max with |y max/min | ∼ ln P 2 ⊥ /q 2 ⊥ . However, in the actual calculation of c n below, it is convenient to integrate over −∞ < y g < ∞. The difference in c n caused by this approximation is power-suppressed. Second, in the soft emission kernel S g , one approximates k J⊥ = 1 2 q ⊥ − P ⊥ ≈ −P ⊥ . Again the difference is powersuppressed in q ⊥ /P ⊥ . There are also power corrections from the hard part that can affect azimuthal asymmetries. In this paper, we do not study these corrections systematically (they are in any case beyond the leading TMD factorization formalism), and in most of our calculations below, we neglect the q ⊥ -dependence of c n . However, in Section IIIA, we will include part of power corrections for phenomenological reasons. When calculating the Fourier coefficients c n , we need to subtract in the k g -integral the configuration where the soft gluon is emitted inside the jet cone of radius R. Namely, we have to impose the constraint As a result, {c n } depend on R rather strongly. To gain analytical insights into the coefficients, it is convenient to replace (7) by which is equivalent to (7) when R 1. We can then obtain the following explicit expression for an arbitrary Fourier coefficient c n c n = ln 1 with The collinear singularity is isolated in the logarithm ln 1/R 2 , and the remaining part is finite. In particular, c 0 = ln(1/R 2 ) and the first few coefficients of the rest read f (1) = 2 ln 4−2, f (2) = −1, f (3) = 2 ln 4−14/3 and f (4) = −5/2. For sufficiently large values of n, we find f (n) ln(b 2 0 /n 2 ) with b 0 = 2e −γ E (γ E is the Euler constant). Also note that g(nR) ≈ n 2 R 2 /4 when nR 1, while g(nR) ≈ ln(n 2 R 2 /b 2 0 ) in the limit nR 1. This indicates that c n vanishes when nR 1. When R is large ∼ O(1), we should return to (7). The Fourier coefficients can be evaluated numerically as follows (see (A3)) where y ± = ± R 2 − φ 2 . For example, for R = 1, we have c 0 −0.25, c 1 = 0.78 and c 2 = −0.30. As shown in Fig. 3, c n decreases approximately as ln 1/R 2 for small n values, while oscillations around zero start to appear for large-n coefficients. We now extend the above one-loop results to all orders in the TMD framework by resumming the double and single logarithms in Q 2 /q 2 ⊥ . This is appropriately carried out in the Fourier transformed b ⊥ -space. The resummed azimuthal averaged cross section reads [33], where Here and in the following, we neglect the high order corrections to the hard factor in the resummation formulas. The Sudakov form factor is defined as To derive the resummation result for the azimuthal angle dependent differential cross section, we first compute the Fourier transform of the soft gluon radiation contribution at one-loop order from Eq. (4), by applying the Jacobi-Anger expansion, and the integration formula, Importantly, the q ⊥ -integral gives a constant although originally in momentum space the angular dependent terms are singular 1/q 2 ⊥ , see, Eq. (4). At higher orders there are double logarithmic corrections but they can be resummed together with the angular-independent term [36,37]. After this resummation, we arrive at An important feature of the above result is that the Fourier coefficients scale as in the small-q ⊥ region [37].
To evaluate (17), following Ref. [92] we employ the socalled b * -prescription to suppress the large-b ⊥ region and introduce non-perturbative form factors associated with the initial and final state radiations, The form factor associated with the incoming quark is [93,94] with Q 2 0 = 2.4 GeV 2 and that for the final state jet is assumed to be We note that there is no constraint for g Λ from experimental data so far. For an illustration, we employ the value g Λ = 0.1 GeV 2 . The numerical results for cos φ, cos 2φ and cos 3φ azimuthal asymmetries in typical kinematics of EIC are presented in Fig. 4 1 . One can see the scaling (18) in the small-q ⊥ region. For narrow jet (R = 0.4, top panel) production, the cos φ modulation is dominant as expected, though the cos(2φ) and cos(3φ) modulations are not negligible. Interestingly, the latter two flip signs for fat jet (R = 1, bottom panel) production, while the cos φ modulation is relatively unaffected. This can be understood from the numerical results of c n (R) as shown in Fig. 3. When the cone size R increases from 0.4 to 1, c 1 remains positive, while c 2 and c 3 become negative.
We close this subsection with a remark on the QED radiative contribution to the azimuthal angle asymmetries. When computing the graphs with a soft/collinear photon emitted from the final state electron, the fixed order calculation produces a large logarithm ln Q 2 m 2 e (roughly ≈ 19 for typical EIC kinematics, see section V), which compensates the smallness of α em to a large extent. This contribution can be considered as part of QED radiative corrections, similar to that discussed in Ref. [97] for inclusive DIS.

B. Photon plus jet production in pp Collisions
Next, we consider photon plus jet production in pp collisions. The dominant partonic channel is q(p 1 )g(p 2 ) → q(k J )γ(k γ ). The leading order cross section of this process is given by, where σ qg→γq −ŝ u −û s with the usual Mandelstame variables for the 2 → 2 partonic processes: In the above equation dΩ = dy J dy γ d 2 P ⊥ d 2 q ⊥ represents the phase space of the final state photon and jet, and y γ and y J are their rapidities. The parton momenta fraction are fixed according to x q,g = P ⊥ (e ±y J + e ±yγ )/ √ s. At one-loop order, the soft gluon radiation gives the follow- 1 For convenience, all cos(nφ) asymmetries in Figs. 4-6 have been multiplied by a factor of n, i.e., the curves in these figures correspond to (n cos(nφ) ).
ing contribution where we used the results in Appendix B. c n are the same as in the previous subsection. Again the singularities in the azimuthally symmetric part can be resummed to all-orders in the TMD framework. Considering that the initial state consists of a quark and a gluon, we obtain the resummed cross section where φ b is the angle between b and P ⊥ . Clearly, the first term in the above contributes to the azimuthal angle averaged differential cross section. The second term contributes to various cos(nφ) asymmetries, which can be further written as, The perturbative Sudakov factor is given by For the numerical evaluation, we need to introduce non-perturbative form factors similarly to (19) and (20). In the present case, we take where Sud q NP is the same as (20) with Q → P ⊥ . The results for cos(nφ) (n = 1, 2, 3) are shown in  Similarly, for the qq → gγ channel with the Born cross section one can obtain the following eikonal factors due to soft gluon emissions which is related to the qg → qγ channel through the crossing symmetry (k J ↔ p 2 ). Thus the corresponding perturbative Sudakov factor reads (30) Similar numeric results can be obtained for this channel as well. However, at the RHIC and LHC kinematics, qq → γg channel is negilible as compared to the qg → γq channel.
A few additional remarks are in order before we leave this section. First, as demonstrated in the above two cases, in which only one of the final state particles is colored while the other measured particle is color-neutral, the coefficients of odd harmonics are non-vanishing, and the dominant azimuthal angle correlation between q ⊥ and P ⊥ is of the form cos(φ). This is simply due to the fact that the soft gluon radiation close to the measured jet is favored. The asymmetry grows with q ⊥ and can easily reach 10% or even > 20% when R is small. Second, we expect that similar conclusions should hold in other processes such as Higgs (or Z/W boson) plus jet production in pp collisions. For the case of the production of Z/W plus jet [34,35,95], the Fourier coefficients c n and the Sudakov coefficients are identical to those in the photon plus jet case. For the case of the Higgs plus jet, the asymmetries are proportional to the same c n , and the angle-independent results including the Sudakov factor can be found in Ref. [96].
At last, since the pure initial-state gluon radiations (represented by S g (p 1 , p 2 ) in our calculation) do not generate asymmetries, and one of the final-state particles is colorless, one can find that only one set of the Fourier coefficients c fi n (see Appendix A and B) arises from the gluon emission related to the final state jet. In contrast, a set of more complicated asymmetries depending on the rapidity difference of the final-state jets can arise from the eikonal factor S g (k 1 , k 2 ), which involves two finalstate particles. In the following section, we consider the production of two colored particles (jets) and study their angular correlations.

III. DIJET IN THE FINAL STATE
In this section we consider diffractive dijet production in γp or γA collisions and inclusive dijet production in γp and pp collisions. Since there are two colored objects in the final state, the dominant asymmetry is expected to be the cos(2φ) term arising from soft gluon radiation along the nearly back-to-back jets. We will only focus on cos(2φ) asymmetry in this section. The extension to other higher harmonics (such as cos(4φ)) should be straightforward.

A. Diffractive Dijet Production
We first study the diffractive photoproduction of dijets, γp → qq + p, γA → qq + A. In this process, an on-shell photon fluctuates into a quark-antiquark pair which then scatters off the nucleon/nucleus target via a color-singlet exchange and forms a final state dijet with momentum k 1 and k 2 . This process can be studied, for example, in ultraperipheral pA and AA collisions at RHIC and the LHC, and the first data from the CMS collaboration came out recently [39]. It can also be studied at the planned EIC in the future. Because the initial state does not carry color, there will be only final state radiation from k 1 and k 2 . Therefore, the soft gluon radiation kernel is simply given by the eikonal factor The above factor corresponds to the classical eikonal radiation from the fast-moving external currents k i , which is valid in the soft limit. Integrating over the phase space (see Appendix B), we can write Since the dijet configuration is symmetric, there is no cos φ term. Note that, for once, we have included q ⊥dependence in the coefficients c diff n for a phenomenological reason (see below). As already mentioned in Section IIA, in general c n depends on q ⊥ due to power corrections. Fourier transforming (32) to the b ⊥ space and including the virtual terms, we find, at one-loop order, After resumming the logarithms and Fourier transforming back, we obtain where We emphasize that the value c diff 0 (q ⊥ = 0) appears in the Sudakov form factor. This is because only the leading power contribution can be resummed into an exponential. To carry out the b ⊥ -integral, we include nonperturbative Sudakov factors. Since there is no TMD quark or gluon distribution involved, we use (see (21)) where the factor of 2 is because there are two jets in the final state. Let us evaluate (34) in two different ways. First we ignore power corrections. In this case, and in the limit a 0 depends on the rapidity difference ∆y 12 = |y 1 − y 2 | as a 0 =ŝ 2 /tû = 2 + 2 cosh(∆y 12 ). The function a 2 (∆y 12 ) remained undetermined in our previous publication [38], but here we can report a fully analytic result ln a 2 = ∆y 12 sinh ∆y 12 − cosh ∆y 12 ln [2 (1 + cosh ∆y 12 obtained after a rather lengthy calculation outlined in Appendix B. Eq. (38) has a very mild dependence on ∆y 12 . It increases slightly from a 2 = 1/4 to a 2 = 1/e when ∆y 12 goes from 0 to ∞ [38]. When R is not very small, we can calculate c diff n numerically using the formula where y + = R 2 − φ 2 as before and we assumed ∆y 12 = 0 for simplicity. In particular, we find for R = 0.4, to c diff 0 , which is a constant in the soft gluon limit. In the small-q ⊥ region, we find cos(2φ) ∝ q 2 ⊥ as expected, and in the large-q ⊥ region the asymmetry reaches a plateau.
However, recent experimental data from the CMS collaboration [39] show a monotonically increasing behavior in the large q ⊥ region. We have noticed that this discrepancy can be alleviated, at least qualitatively, by including the following two sources of power corrections. First, when q ⊥ gets larger, |k 1,2⊥ | differ from |P ⊥ | = |k 1⊥ − k 2⊥ |/2 significantly. We can correct for this difference by computing the soft emission kernel using the exact kinematics where k 1,2⊥ = q ⊥ /2± P ⊥ , φ 1 , φ 2 and φ g are the azimuthal angles of k 1⊥ , k 2⊥ and k g⊥ with respect to P ⊥ , respec- tively. Second, we impose precise rapidity cutoffs for the y g -integral. For instance, when These power corrections effectively make c diff 0 and c diff 2 dependent on q ⊥ . It is interesting to note that c diff 0 is mainly affected by the first source, whereas c diff 2 receives contributions from both sources with opposite signs. We have re-evaluated (34) taking this q ⊥dependence into account. The result is that cos(2φ) now becomes an increasing function in the large-q ⊥ region as shown in the lower panel of Fig. 7. The wiggles on the leading order curve are caused by the cancellation of the two sources in c diff 2 mentioned above. Of course, this is not a fully consistent procedure as we ignore other possible sources of power corrections such as those coming from the hard part. Yet, the better agreement with the CMS data may suggest that these are an important part of power corrections. In Fig. 8, we plot the asymmetry for EIC kinematics. Its size and q ⊥ dependent behavior is similar to that in the CMS kinemic region except that the non-perturbative effect is more significant.
We can also derive the resummed, full φ-dependent cross section. Instead of Fourier expanding as in (32), we can decompose the soft factor as where the term in the bracket is free from infrared divergence after averaging over φ. The last term can be combined with the virtual contribution and get exponentiated. The resummed soft factor thus takes form, where f (b ⊥ ) is the Fourier transform of To evaluate (43) efficiently, we avoid the task of computing the Fourier transform of f (k ⊥ ). Instead, we proceed by rewriting (43) as where The subscripts s and a denote the azimuthally symmetric and asymmetric parts, respectively. The k ⊥ -integral in (45) can then be done straightforwardly after including the nonperturbative Sudakov factor.
We plot the resummed azimuthal angle distribution of q 2 ⊥ S(q ⊥ ) in Fig. 9 (lower panel) for different values of q ⊥ , and also the leading order result (upper panel) for comparison . The scale of the coupling constant in the leading order calculation is chosen to be P ⊥ . One can see that the φ distribution becomes smoother after performing the all-order resummation.
B. Inclusive Dijet in γp Collisions through γg → qq Next we turn to inclusive dijet photoproduction γp → jjX. Here we focus on the direct photon contribution through γg → qq channel. The leading order cross section of this process is given by, where σ γg 0 represents the leading order cross section, dΩ = dy 1 dy 2 d 2 P ⊥ d 2 q ⊥ for the phase space, and f g (x g ) is the gluon distribution with x g = P ⊥ (e y1 + e y2 ) / √ s γp is momentum fraction of the nucleon carried by the gluon. The amplitude squared of the one soft gluon radiation can be written as [28], where M 0 represents the leading Born amplitude and the expression inside the square brackets represents the colorweighted sum of the eikonal radiation functions. All pairs external color lines need to be summed over. Applying the results from the Appendix, we obtain the soft gluon radiation contribution to the differential cross section, In the small-R limit, we have where a 1 = 1/e and a 2 is the same as in the previous subsection. In the following numeric calculations, we will apply c γg 0 = 3.14 and c γg 2 = 0.96 for R = 0.4 and the two jets are at the same rapidity. We now perform the resummation of double logarithms in the standard TMD framework, where Again, we separate the Sudakov form factor into the perturbative and non-perturbative parts The perturbative part is given by and we apply the b * -prescription. The non-perturbative part for the present problem is with Q → P ⊥ in (20). In Fig. 10, we show the numerical results for cos(2φ) as a function of q ⊥ for the typical kinematics at the future EIC, √ s γp = 100GeV, P ⊥ ∼ 15 GeV and the two jets are at the same rapidity. Compared to the results in the diffractive case studied in the above subsection, the impact of resummation is more pronounced.
C. Inclusive Dijet in pp Collisions from gg → gg channel At the LHC, dijet production is dominated by the gg → gg channel. The soft gluon radiation contribu-tion to the azimuthally symmetric part of the differential cross section has been derived in Ref. [30]. In this subsection, we extend this work to the angular dependent part. The leading Born amplitude can be decomposed as where a, b, c, d are color indices in the reaction p 1a p 2b → k 1c k 2d . From the above color structure, we notice that A 1,2,3 represent gauge invariant amplitudes in the s, t, uchannels, respectively. After summing over gluon helicities, one finds the following useful results where N represents the overall normalization. The leading order amplitude squared can be written as which is consistent with the well-known gg → gg Born amplitude square. The soft gluon radiation amplitude soft at one-loop order (real diagram) takes the form [30], Squaring this, we get By applying the following relation, we can rewrite the amplitude squared as, Furthermore, we use the following results for A 1 , A 2 and A 3 , We emphasize that the above results are gauge invariant components in the amplitude squared of the gg → gg channel. Using the results of Appendix B, we obtain the differential cross section where σ gg 0 represents the leading order cross section, dΩ = dy 1 dy 2 d 2 P ⊥ d 2 q ⊥ for the phase space, and f g (x g ) is the gluon distribution with x 1 = P ⊥ (e y1 + e y2 ) / √ s and x 2 = P ⊥ (e −y1 + e −y2 ) / √ s are momentum fractions of the incoming hadrons carried by the gluons. In the small-R limit, we have where a 0,1,2 are the same as in the previous sections. The all-order resummation for dijet production in pp collisions at the next-to-leading logarithmic level has to be done in a matrix form in color space [30]. This is because the final state jets and incoming partons form a color antenna with various representations of the color SU(3) group [98,99]. For simplicity, we work in the improved leading logarithmic approximation (LLA ) where we include only the diagonal part in color space, namely, the leading double logarithms and those single logarithms associated with the initial parton distributions and final state jets. (The terms which depend on kinematic variables in c gg 0 are omitted.) In this approximation, we can write Note that only the jet-size dependent term in c gg 0 of Eq. (66) was included in the above Sudakov form factor. The rest should be included in the matrix form of the resummation beyond the approximation we adopted here. For the cos(2φ) term, we have The numerical estimate of the resulting cos(2φ) asymmetry for the LHC kinematics (after including nonperturbative Sudakov factors) has been presented in Ref. [38]. Going beyond the LLA , we need to use the matrix form of resummation. The matrix for the cos(2φ) term may be different from that for the azimuthally symmetric term [37]. For completeness, below we present the results for the other partonic channels in pp collisions. Following the same procedure as shown above for the gluon channel and employing the results summarized in Appendix B, one can obtain the corresponding expression for the qq → qq channel as follows In the small-R limit, one finds where c qq 0 agrees with the result in Ref. [30] with N c = 3. Similarly, the results for the gg → qq can be cast into In the small-R limit, one finds the expressions for the coefficients For the inverse process of the above channel, i.e., qq → gg, one simply can replace the color factors of ln P 2 ⊥ q 2 ⊥ and Fourier coefficients by C F and C A in Eq. (72). Then the first two coefficients are For the qg → qg channel, we have the following soft gluon radiation contribution to the differential cross section, Taking the small-R limit, one finds where c qg 0 also agrees with the result in Ref. [30] with N c = 3 and c qg 2 is new. Here we have neglected the odd harmonics since those terms vanish when final state jets are symmetrized.
Let us comment on the patterns of the above logarithms that one observes from the various processes. First, the color factors of the term ln are associated with the incoming partons, while the color factor of the logarithm ln a0,1 R 2 are determined by the final state jets. The second term of the coefficients c 0 and c 2 are process dependent. In particular, the ln a2 a1 terms are expected to be small and it vanishes when ∆y 12 → ∞.

IV. DIJET PRODUCTION IN DIS TO PROBE THE LINEARLY POLARIZED GLUON DISTRIBUTION
In this section, we return to inclusive dijet production in DIS off a nucleon/nucleus Differently from Section IIIB, here the exchanged photon is virtual with invariant mass squared q 2 = −Q 2 . As mentioned in the Introduction, this process has been proposed to study one particular aspect of the gluon distribution in the nucleon/nucleus, the so-called linearly polarized gluon distribution [17][18][19][20]100]. The dependence of the differential cross section on q ⊥ = k 1⊥ + k 2⊥ is sensitive to the TMD gluon distributions of the nucleon/nucleus. Among them, the linearly polarized gluon distribution will lead to a characteristic cos(2φ) asymmetry, where φ is the azimuthal angle between q ⊥ and P ⊥ [18]. This observation has gained more importance when it was realized that the linearly polarized gluon distribution is of the same size as the usual gluon TMD distribution in the small-x saturation formalism [19]. Since then, several proposals have been made to measure the cos(2φ) asymmetry at the planned electron-ion collider [17,20]. Moreover, the distribution has been widely applied to many other processes [21,22,[101][102][103][104][105][106][107]. However, there are two important issues which complicate the interpretation of such measurements, but have not been adequately investigated in the literature. First, the collinear gluon radiation from the incoming parton can generate the cos(2φ) modulation. This has been well understood in the collinear factorization framework at moderate transverse momentum and also in the TMD resummation formalism [48,49,92,108,109]. This perturbative effect can mimic the nonperturbative, intrinsic modulation due to the linearly polarized gluon distribution, but it has been largely ignored in the previous phenomenological studies. Another source of the cos(2φ) correlation that has been missing in the literature of the linearly polarized gluon distribution is the soft gluon radiation from the final state jets as pointed out in [38] and discussed in the previous sections. Therefore, in order to reliably extract the linearly polarized gluon distribution through the measurement of the cos(2φ) asymmetry, it is important to quantify these 'background' effects. 2 In this section, we perform a systematic study of the cos(2φ) asymmetry in DIS dijet production including the above three physics: the 'intrinsic' and 'collinear radiation generated' linearly polarized gluon distributions, and the final state soft gluon radiation contribution. We shall focus on the TMD domain, i.e., P ⊥ q ⊥ , where the leading jet transverse momentum is much larger than the total transverse momentum of the two jets. In this region, we have to perform an all-order resummation of the logarithms (α s ln 2 P 2 ⊥ /q 2 ⊥ ) n . Previously, the TMD resummation for the linearly polarized gluon distribution has been studied in Refs. [48,49,108] and its impact on Higgs Boson production was found to be very small [112,113]. Resummation effects on the cos(2φ) in photon-jet correlation have also been studied in Ref. [21], which, however, only included the intrinsic linearly polarized gluon contribution. Our calculations in the following will show that the collinear radiation generated linearly polarized gluon distribution dominates over the intrinsic one at higher hard momentum scales Q 2 , say, at the scale of the Higgs mass. However, their relative importance strongly depends on Q 2 , and this leaves an opportunity to explore the transition from intrinsic to collinear radiation regimes in future experiments.

A. Linearly Polarized Gluon Distribution
In this subsection, we briefly introduce the linearly polarized gluon distribution and study the associated QCD evolution. The TMD gluon distributions are defined through the following matrix element [114][115][116], where the nucleon moves along +ẑ-direction and F µν a is the gluon field strength tensor. The light-cone components are defined as k ± = (k 0 ± k 3 )/ √ 2. In the above equation, x is the longitudinal momentum fraction carried by the gluon and k ⊥ is the transverse momentum. The gauge link L v is constructed in the adjoint representation and depends on the process [117]. For an unpolarized nucleon at leading twist, the above matrix element contains two independent TMD gluon distributions [115], where g µν ⊥ has only transverse components g ij ⊥ = δ ij . f g (x, k ⊥ ) is the usual azimuthally symmetric TMD gluon distribution, and h g (x, k ⊥ ) is the linearly polarized gluon distribution. h g vanishes when M µν is integrated over transverse momentum, which means there is no integrated version of the linearly polarized gluon distribution h g (x).
To apply TMD gluon distributions in hard scattering processes, we have to take into account the TMD evolution and resummation. For the azimuthally symmetric part, we follow the "standard" scheme (also called Collins 2011 Scheme) [94,[118][119][120] which reads, in coordinate space, where the perturbative Sudakov factor is processdependent (to be specified below) and the nonperturbative part is given by Sud g We have set the rapidity regulator ζ c and the renormalization scale µ 2 to be both Q 2 , and F g (α s (Q)) = 1 + O(α s ) in the standard TMD scheme. For the C-coefficients, we use the one-loop results in the standard TMD scheme. Note that in this scheme C g/g vanishes at one-loop order. Numerically, we find that the contribution from C g/q is negligible. We therefore only keep the delta function term in (82) in the following.
The Collins-Soper evolution equation for the linearly polarized gluon distribution can be derived in a similar manner. Again in the b ⊥ -space, we parametrize as where the tensor structure is uniquely determined by the traceless condition. The solution of the evolution equation takes the form The distribution at the lower scale ζ c = µ b does not contain large logarithms, but it cannot be written similarly to (81) because there is no integrated h g distribution.
In the large-b ⊥ region where physics becomes nonperturbative, it has to be modeled. Following Ref. [102], we parameterize the linearly polarized gluon distribution in term of the normal gluon distribution, where Q h ≈ 1 GeV and f g (x, µ) is the integrated gluon distribution. In this model, we assume that the linearly polarized gluon distribution has the same x-dependence as the normal gluon distribution f g (x, µ). In reality, they may be totally different. Eq. (85) with (86) (extrapolated to the full b ⊥ region including small-b ⊥ ) is what we call the intrinsic part of the linearly polarized gluon distribution. In momentum space, it is proportional to k 2 ⊥ in the small k ⊥ region and satisfies the positivity bound h g < f g . As mentioned already, in most literature only this part has been used to calculate the cos(2φ) asymmetry. In such approaches,h g (b ⊥ ) ∼ b 2 ⊥ as b ⊥ → 0. However, we know that the small-b ⊥ behavior of the linearly polarized gluon distribution is perturbatively calculable via collinear gluon radiation at large transverse momentum The C-coefficients start at O(α s ), and are given by [48,49,108,109] The two-loop results have been recently derived in Ref. [109]. In the following, as an illustration we only use the one-loop results. We see from (87) that the correct small-b ⊥ behavior is not quadratic but constant (up to the logarithmic running of the coupling) This is so because there is no dimensionful parameter in perturbation theory to compensate the dimension of b 2 ⊥ . Comparing the large and small b ⊥ behaviors of the intrinsic (86) and radiative (87) contributions, we notice that they dominate h g (x, b ⊥ ) separately in these two regions. Therefore, we can combine them together and arrive at the following two-component model At large b ⊥ , the first term dominates, while at small-b ⊥ , the second term dominates and the first term is negligible. An important feature of (91) is that the scale (Q 2 ) dependence of the linearly polarized gluon distribution is dictated by the Sudakov form factor. At very large Q 2 , it pushes the b ⊥ -distribution to the small-b ⊥ region, and we only need to take into account the contribution from the second term. At low Q 2 , the large-b ⊥ region will be important to obtain the k ⊥ distribution, and one may keep only the first term to a good approximation. For moderate values of Q 2 , we may need to take into account both contributions. Therefore, by studying different hard processes sensitive to the linearly polarized gluon distribution, we will be able to investigate the transition from non-perturbative to perturbative regimes. This provides a unique perspective for the nucleon/nucleus tomography study in future experiments, in particular, at the planned EIC.
In Fig. 11, we show the ratio h g (x, k ⊥ )/f g (x, k ⊥ ) as a function of k ⊥ for different values of Q 2 at a fixed value of x = 0.05. These plots clearly demonstrate the above points. In particular, it is interesting to notice that, for Higgs boson production at the scale Q = M h , the linearly polarized distribution is completely dominated by the collinear radiation contribution. On the other hand, when Q = 5 GeV, the k ⊥ -dependence is dominated by the non-perturbative part. Of course, we do not know the actual magnitude of the nonperturbative part. (Eq. (86) is just a model.) However, our results suggest that this can be constrained by scanning Q 2 in the relatively low momentum region at the EIC.
In the above discussions, the separation of the "intrinsic" and "collinear" parts of the linearly polarized gluon distribution depends on the model we used, where they have different behaviors at large and small b ⊥ . It will be interesting to develop a model to capture both features in a single set-up. In particular, in the small-x dipole formalism, both the linearly polarized gluon distribution and the normal gluon distribution can be calculated from the same dipole amplitude [19] and it may be possible to include the "intrinsic" and "collinear" parts at the same time. Further developments are needed to implement collinear gluon radiation contribution in the gluon distribution functions in the small-x dipole formalism, see, for example, the discussions in Ref. [121].

B. cos 2φ Correlation in Inclusive Dijet Production in DIS
We now calculate the cos(2φ) asymmetry in inclusive dijet production in DIS using the linearly polarized gluon distribution constructed in the previous subsection. We focus on the gluonic channel γ * T,L g → qq where the incoming photon can be transversely (T ) or longitudinally (L) polarized. In the collinear factorization framework, the differential cross section can be written as, where dΩ = dy 1 dy 2 d 2 P ⊥ d 2 q ⊥ . f g represents the integrated gluon distribution. At the leading Born order, we haveσ where ξ = x g /x with x g being the momentum fraction carried by the gluon. It can be determined from the dijet kinematics as x g = ( where s is the center of mass energy of the γ * p system and z is the momentum fraction of the virtual photon carried by the quark jet. Gluon radiations from the incoming gluon and the outgoing qq pair will generate not only a nonzero transverse momentum q ⊥ but also a cos(2φ) asymmetry. In the TMD kinematics P ⊥ q ⊥ , we have both collinear and soft gluon contributions, where σ 0 and σ 2 are normalization factors for the differential cross sections and c γg 0 and c γg 2 are the same as those defined for γg → qq subprocess in Sec. IIIB. P (<) g/g denotes the collinear splitting kernel without the delta function part. The soft radiation part is essentially the same as (50) except that now the phase space has increasedŝ →ŝ + Q 2 since the incoming photon is virtual. Note that the Mandelstam variables are related bŷ stû = P 2 ⊥ (ŝ + Q 2 ) 2 in this case. The above equations apply to both transverse and longitudinal incoming photons. To leading order, we have the following relations for the cross section ratios [19,100], The singularity q ⊥ → 0 can be factorized and resummed in the b ⊥ -space into the TMD gluon distributions f g and h g , as well as the soft factors associated with the final state jets. This converts (92) into where the azimuthal symmetric term can be written as The perturbative Sudakov form factor is defined as and the nonperturbative part is the same as (56).
On the other hand, the cos(2φ) term consists of two parts The first term comes from the soft gluon emission from the final state jets [38], and the second term comes from the linearly polarized gluon distribution whose resummation has been discussed in the previous subsection. 3 Note that the latter is absent in photo-production studied in [38] and also in Section III because σ 2 vanishes when Q 2 = 0, see (96). All in all, we arrive at the following representation Eq. (101) clearly exhibits the three distinct contributions 3 Our one-loop result only demonstrates the contribution from the linearly polarized gluon distribution. The associated Collins-Soper evolution for the linearly polarized gluon distribution will eventually lead to a resummation following the discussion in previous section. There is also a soft factor associated with the final state jet. We expect the same resummation formalism as that in W 0 . mentioned at the beginning of this section: the intrinsic (∼ g h ) and radiative (∼ C h/i ) contributions to the linearly polarized gluon distribution, and the soft gluon emission contribution (∼ c γg 2 ). The numerical results for the cos 2φ asymmetry are presented in Fig. 12 for longitudinal (top panel) and transverse (bottom panel) virtual photons for a typical EIC kinematics with R = 0.4. While the contribution from the linearly polarized gluon distribution (dashed curve) is noticeable, it is overwhelmed by that from the final state soft gluon emissions in the whole range of q ⊥ . The latter is independent of the polarization (longitudinal/transverse) of the virtual photon. In order to extract h g from this observable, it is probably better to use larger values of R, say, R = 1 to suppress the final state emissions. If one is ultimately interested in the 'intrinsic' part of h g , further considerations are required (such as lowering Q) to suppress the 'collinear radiative' part of h g , see Fig. 11.

V. LEPTON PAIR PRODUCTION IN TWO-PHOTON PROCESS
As a final example, we consider lepton pair production in QED γγ → + (k 1 ) − (k 2 ) which has been actively studied recently in ultraperipheral collisions (UPC) at RHIC and the LHC [66][67][68][69][70][71][72]. Similarly to the dijet problem in QCD, the dilepton azimuthal correlation is dominated by soft photon radiations from the final state leptons at small q ⊥ = k 1 + k 2 , where M (0) represents the leading order Born amplitude and the soft photon carries momentum k s . Working in the laboratory frame, we integrate the soft emission kernel over the photon rapiditŷ where m 2 t ≡ m 2 /P 2 ⊥ and m is the lepton mass. Because of the mass term, there is no collinear divergence associated with final state radiations. In other words, m t plays the role of R in the previous sections. The integral is carried out in Appendix C, with the following result for the oneloop soft factor where Q 2 the invariant mass squared of the lepton pair.
Adding the virtual contribution, we obtain the soft factor Comments are in order regarding the relation to Ref. [81]. In this reference, the photon rapidity integral was carried out in the lepton frame (outgoing leptons are along the z-axis), not in the lab frame. This resulted in a different expression than (104), and the two results lead to slightly different predictions for the acoplanarity of the lepton pair in UPC at the LHC, mainly in the moderate acoplanarity region. However, both predictions are compatible with the data due to large experimental uncertainties.
In addition to the angular independent part (104), we have also calculated the cos(2φ) term in Appendix C. In the small mass limit, the ratio c 2 /c 0 is close to unity. This is because soft photons are concentrated around the lepton directions due to the collinear enhancement. In the limit m → 0, the φ-distribution of photons diverges around φ = 0 and π.
All-order resummation can be carried out similarly to the previous problems. In place of the nonperturbative Sudakov form factor, we introduce a simple Gaussian factor which takes into account the intrinsic transverse momentum of incoming photons p ⊥ . The total transverse momentum distribution can be written as where β = αe π ln Q 2 m 2 . For simplicity, let us set y 1 = y 2 and define The general results provided in Appendix C can be used to compute cases with arbitrary ∆y. We employ the value Q 0 = 40 MeV ∼ 1/R A for the Gaussian width. The Fourier transform leads to the following result where 1 F 1 is a hypergeometric function. We thus arrive at For the typical kinematics of dimuon production from the two photon processes in UPC heavy ion collision at the LHC, β ≈ 0.02 and m/P ⊥ ≈ 0.02 or c 2γ ≈ 0.70. For RHIC kinematics, the typical values of P ⊥ range from 350 MeV to 2 GeV. Taking, for example, Q = 2P ⊥ = 1 GeV (y 1 = y 2 ), we get c 2γ ≈ 0.82 and β ≈ 0.035 for di-electron production. In Fig. 13, we plot cos(2φ) as a function of q ⊥ /Q 0 for the above RHIC and LHC kinematics. We find that the asymmetry is very small when q ⊥ is below 2Q 0 . However, it rapidly grows around 2Q 0 and becomes sizable at 3Q 0 which is about 100 MeV. Similar results can be obtained for cos(4φ) as well. Both the cos 2φ and cos 4φ azimuthal asymmetries in di-lepton production can also be induced by the primordial linearly polarized photon distribution [69,76,77]. Our result shows that the high-q ⊥ tail of the asymmetries are overwhelmingly developed via perturbative final state soft photon radiations. In contrast, at low q ⊥ the asymmetries are mainly attributed to the primordial linearly polarized photon distribution.

VI. CONCLUSIONS
In summary, we have performed a systematic study of azimuthal angular correlations between the total q ⊥ = k 1⊥ + k 2⊥ and relative P ⊥ = ( k 1⊥ − k 2⊥ )/2 transverse momenta in dijet and related systems. The resummation of the double and single logarithms (ln P ⊥ /q ⊥ ) n developed in [36][37][38] is demonstrated for a number of processes. The dominant Fourier modes are cos(2φ) and cos(φ) for dijet and single jet (plus a color-neutral particle) productions, respectively. The expectation value cos(nφ) grows as q n ⊥ in the small-q ⊥ region, and it can easily reach 10-20% in the large-q ⊥ region. While this is an interesting feature of soft gluon emissions in its own right, it can become a serious background for certain purposes. In particular, the cos(2φ) asymmetry in inclusive dijet production proposed as a signal of the linearly polarized gluon distribution [18][19][20][21][22] inevitably faces this challenge.
A majority of the processes we have studied are relevant to the future EIC experiment [9][10][11], such as leptonjet correlation in DIS (Sec. IIA), diffractive (Sec. IIIA) and inclusive (Sec. IIIB) dijet production in photonproton collisions, and inclusive dijet production in DIS (Sec. IV), see the plots in Figs. 4,8,10,12. The comparative study of all these processes will provide a crucial test of our predictions and lead to a better constraint on the nucleon/nucleus tomography in terms of the gluon Wigner distribution and the linearly polarized gluon distribution.
One of the important directions for future research is the understanding of power corrections of the form (q ⊥ /P ⊥ ) n . We have heuristically noticed in Section IIIA (see Fig. 7) that the inclusion of power corrections in the soft emission kernel results in a better agreement with the CMS data in the large-q ⊥ region. The importance of power corrections in the hard part for the extraction of the linearly polarized gluon distribution has been discussed recently [110,111]. A combined analysis of resummation and power-corrections seems to be necessary to correctly interpret the experimental data. In this appendix, we evaluate the integral (4) relevant to lepton-jet production. Introducing the rapidities y g , y J of the soft gluon and the ougoing jet, respectively, we can write where ∆y ≡ y g − y J . Λ is the rapidity cutoff which can be determined by the kinematical constraint k + g < p + 1 , or more explicitly, k g⊥ e yg < p + 1 → e ∆y < p + 1 e −y J k g⊥ → ∆y < Λ ≡ 1 2 ln where Q 2 = −t = −(k 1 − p 1 ) 2 = p + 1 k J⊥ e −y J . Here we define p ± ≡ p 0 ± p 3 and parametrize the momentum of a particle as p µ = (p + , p − , p ⊥ ) where p µ p µ = p + p − − p 2 ⊥ . For massless particles, p ± = p ⊥ e ±y . After integrating over ∆y assuming Λ 1, one finds where only the leading power contributions are kept. We can also implement the constraint (7) by restricting the y g integral for R > |φ| as where y ± = ± R 2 − φ 2 . The result can then be expanded in Fourier series S g (k 1 , p 1 ) ⇒ 1 q 2 ⊥ ln Q 2 q 2 ⊥ + ln Q 2 k 2 J⊥ + c 0 + 2c 1 cos(φ) + 2c 2 cos(2φ) + · · · . (A5) 2 → 2 processes a(p 1 ) + b(p 2 ) → c(k 1 ) + d(k 2 ) with four on-shell external massless particles as follows c fi n cos(nφ)(−1) n , (B11) c fi n cos(nφ)(−1) n , (B13) where we find in the small cone limit c fi 0 = ln 1/R 2 , c fi n = ln 1/R 2 + f (n), c ff 0 = ln 1/R 2 , and It is clear that S g (p 1 , p 2 ) is independent of azimuthal angle since the pure initial state gluon radiation is expected to be symmetric. In general, there are two types of anisotropy generated from the final state radiations, as shown above. The eikonal factors S g (k i , p j ) involve one final-state jet and one initial-state particle, and their contributions are captured by the coefficients c fi n . For odd Fourier coefficients, there is a sign change between S g (k 1 , p i ) and S g (k 2 , p i ). This is due to the fact that the final state gluon radiation is favored along the jet direction and it contributes to odd coefficients oppositely (as shown in Fig. 1 cos(π − φ) = − cos φ) for two jets that are back-to-back. Besides, more complicated angular correlations characterized by c ff n can arise from S g (k 1 , k 2 ) that depends on both final-state jets. Note that in e + q → e + jet scattering in DIS,ŝt/û = Q 4 /k 2 J⊥ so (B14) reduces to (A5) when p 2 and k 2 in Eq. (B14) are identified as the incoming and outgoing quark momenta, respectively. Together with the usual factor of α s /(2π 2 ) and the corresponding color factor, the above identities give rise to the one-loop results for various scattering processes.