The suppressions of dijet azimuthal correlations in the future EIC

Quark-antiquark pair (or dijet) production at the electron-ion collider (EIC) has been argued to be one of most important processes that allowing to access the Weizs\"acker-Williams (WW) gluon distributions at small $x$ limit. Within the framework of Color Glass Condensate (CGC) effective field theory (EFT), we calculated the dijet cross sections and the azimuthal correlations by including the Sudakov resummations, numerical results shown that the back-to-back correlations are significantly suppressed when the Sudakov resummations are taken into account. In addition, by using the solutions of running-coupling Balitsky-Kovchegov (rcBK) equation, the unpolarized and linearly polarized WW gluon distributions both in coordinate and momentum space are presented.


I. INTRODUCTION
Exploring the multi-dimensional structures and the detailed dynamics of a hadron has been one of the most primary goals in hadronic physics. In recent years, the so-called transverse momentum distributions (TMDs) associated with their evolutions, representing essential aspects of the partonic structures in momentum space (parton distributions not only as a function of longitudinal momentum fraction x, but also involves the intrinsic transverse momentum k ⊥ ), have received a lot of concerns both in theory and experiment sides. A future planed deep inelastic scattering (DIS) facilities, such as electron-ion colliders (EIC) [1][2][3], LHeC [4,5] and EicC [6], will provide the cleanest environment to determine these fundamental quantities over a broad range of kinematics at high accuracy.
The quark-antiquark pair (or dijet) production at DIS has been argued to be an important process that allowing to access the partonic structures in momentum space. Especially, when the produced jets are mostly back-to-back, the novel process can be used to as a proxy to probe the Weizsäcker-Williams (WW) gluon distribution [7]. The WW gluon distribution, as one type of gluon TMDs, which features an intrinsic gauge link structures: either future-pointing or past-pointing gauge links in light-cone gauge and can be interpreted as the number density of gluons inside the scattered target [7][8][9]. As addressed in Ref. [10], a nontrivial feature has been realized that these gluons can be linearly polarized even if the hadron is unpolarized, especially, if sizable, can give a large impact on physical observables, for examples, the (pseudo)scalar particle productions in hadronic collisions [11][12][13], as well as the azimuthal correlations of dijets in DIS [7,[14][15][16][17][18][19][20][21][22]. In Refs. [23,24], for dijets production in DIS, it has been shown that unpolarized WW gluon distribution associated with its partner, linearly polarized one, enter the cross section through an azimuthal modulation. Recent developments also shown that the two gluon distributions not only exist in DIS processes, but also wildly appear in other processes, such as the dijets or quark-antiquark pairs [25][26][27], photon pairs [28], muon pairs [16], quarkonium [29][30][31][32] productions in hadronic collisions. A particular feature is that the unpolarized WW gluon distributions are always accompanied by its linearly polarized parts through an azimuthal modulation, this suggests that these gluon distributions (especially for the linearly polarized parts) can be extracted by measuring the azimuthal observables. Another interesting feature has been pointed out that the linearly polarized WW gluon distribution will tend to its saturated bound at large transverse momenta, while is suppressed at small transverse momenta region [23,24], which was then further confirmed by numerical calculations in Refs. [21,22], for recent review please see Ref. [33] and therein. Up to now, the two gluon distributions, apart from the positive bound of linearly polarized part, are still little known experimentally, therefore, to precisely extract and quantitatively understand the behaviors of the two gluon distributions at high energy limit, which allows us in depth understanding the transverse tomography of proton and nucleus, as well as the underlying dynamics at high energy scatterings.
In this paper, we continue to concentrate on the quark-antiquark pair (or dijet) production in a special process that a virtual photon scattering off a hadronic target, this process typically corresponds to the quark-antiquark (or dijet) production in DIS. Due to a broad range of kinematics that can be accessed in the future EIC, this process allows to mapping out the evolutions of gluon TMDs. Especially, when the scattered energy is sufficiently large (or equivalent to small-x), a high density regime can be probed, in this sense, gluon saturation effects will become more predominant, and hence the behaviors of gluon TMDs in the dense regime can be investigated as well. In the correlation limit, prior works have shown that this process can be consistently described in the framework of Color Glass Condensate (CGC) theory, and simultaneously, the relevant two WW gluon distributions can be expressed in terms of Wilson line correlators that making them can be calculated by using the solutions of Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) equation [21][22][23][24]. However, the special limit has a strong restriction on the dijets, which requires the transverse momentum imbalance q ⊥ is much smaller than the relative transverse momentum P ⊥ , in this case, there is another class of effects, known as Sudakov resummations, can give sizable contribution to the cross section, and this should be taken into account when exploring the behaviors of the two gluon distributions in the considered process. Other phenomenological studies have shown that such resummations is very important in quantitatively describing experimental observables [34][35][36][37][38][39][40][41][42][43][44][45][46]. Therefore, in the present work, we plan to conduct this direction by including the Sudakov resummations, and show how the resummations affect the azimuthal correlations of the dijets. Furthermore, as shown in Refs. [21,22,27,47], thanks to the numerical implementation of JIMWLK equation, the two gluon distributions appeared in the considered process can be obtained. However, the numerical implementation is only achieved by fixing coupling constant. The running-coupling corrections, as an important QCD dynamics, has not been performed yet in their studies. In the present work, we prefer to work within the Gaussian approximation of JIMWLK evolution and use the running-coupling Balitsky-Kovchegov (rcBK) equation to take the relevant calculations, because we expect the running-coupling corrections here to be much more important than that of Gaussian approximation. Besides that, we stress the rcBK solutions are well constrained from DIS data [48] and have been received a great successful in describing particle spectrum at hadronic collisions [49,50], therefore, one can be firmly convinced that the gluon TMDs can be well studied by using the solutions of rcBK equation.
The paper is organized as follows. In next section, we will begin with a brief review of the CGC results for the quark-antiquark pair production in DIS, and show that the cross sections are the sensitive probes of gluon TMDs (unpolarized and linearly polarized WW gluon distributions) inside the scattered target at small-x. In Sec. III, we discuss the detailed behaviors of unpolarized and linearly polarized WW gluon distributions in the saturation model and then present the numerical results by using the solutions of rcBK equation. In Sec. IV, we give our numerical results for the dijet cross sections, as well as discuss how the effects of Sudakov resummations and gluon saturation affect their azimuthal correlations. Finally, a short summary is given in Sec. V.

II. DIS DIJET PRODUCTION IN COLOR GLASS CONDENSATE FORMALISM
In this section, we firstly give a brief review of the quark-antiquark pair (or dijet) differential cross section in the CGC formalism, and show that, by applying the "back-to-back correlation limit", the cross sections are expressed in terms of WW gluon distributions.
To begin with, we consider a high energetic virtual photon γ * with four-momentum p scattering off a hadronic target, after that, a quark-antiquark pair qq is produced. This process typically corresponds to quark-antiquark pair (or dijet 1 ) production in DIS via the exchange of a virtual photon. For this process, the hadronic target can be regarded as classical gluon field due to the dominance of gluon saturation effects [52][53][54][55][56]. Therefore, it is convenient to calculate the cross section in the CGC formalism, the standard way is to consider the virtual photon γ * fluctuates into a quark-antiquark pair qq, and then traverse the background field in an eikonal way, simultaneously, the fluctuated quarks receive a contribution to the scattering amplitude with a Wilson line. After averaging over the photon's polarization and summing over the quark's helicities and colors, the differential cross section for this process can be expressed as follow [7,23], where α em is the electromagnetic structure constant, N c and e q are the number of colors and the electrical charge of the produced quark. The two-dimensional variables x 1 and x 2 denote the transverse coordinate of the outgoing quark and antiquark in the amplitude respectively. Clearly, the variable u = x 1 − x 2 and v = zx 1 + (1 − z)x 2 are the transverse size of the produced quark-antiquark pair and the transverse position of the incoming virtual photon, where z = k + 1 /p + is the longitudinal momentum fraction of the outgoing quark with respect to the incoming virtual photon. Similar notations hold for the complex conjugate amplitude. The splitting function ψ T,Lλ αβ (p + , z, u) denotes the virtual photon γ * with longitudinal momentum p + and virtuality Q 2 fluctuates into a quark-antiquark pair with size u, this function can be calculated in light-cone perturbative theory and be expressed in the following form [57], where α, β are the helicities of produced quark and antiquark, λ is the photon's polarization, T, L denote the transverse and longitudinal components of the wavefunction, 2 f = z(1 − z)Q 2 + m 2 q with m q is the mass of produced quark. It should be noted that the interactions between the produced quark-antiquark pair and the scattered target are fully characterized by Wilson line correlators, as shown in the last line of Eq. (1). Specifically, the term equal to 1 means the virtual photon freely pass through the target, and the quark-antiquark pair without scattering in the amplitude and complex conjugate amplitude. The interference terms S 2 xg (x 1 , x 2 ) and S (2) xg (x 2 , x 1 ), known as dipole, denote the quark-antiquark pair scatters off the target either in the amplitude or the complex conjugate amplitude. The four-point function S (4) xg (x 1 , x 2 , x 1 , x 2 ), named as quadrupole, describes the quark-antiquark pair undergoes multi-scattering with target both in amplitude and complex conjugate amplitude. The dipole and quadrupole are defined as: with The notation · · · xg used here denotes the averaging over all possible distributions of color sources and x g being the fraction of longitudinal momentum in the scattered target. As shown in the above, the quark-antiquark pair (or dijet) production in DIS process can be quite well described in the CGC formalism, and the cross section directly probes the internal dynamics of the scattered target. It should be reminded that the CGC approach does not dependent on momentum scales and, in principle, can be applied in any scattering process when high density or "saturation" regime accessed. This treatment is different from the TMD factorization approach, since the later one is strongly dependent on momentum scales. In order to recover the TMD results, it is useful to change variables as follow, then, Eq. (1) casts into, where q ⊥ = k 1⊥ + k 2⊥ and P ⊥ = (1 − z)k 1⊥ − zk 2⊥ are the transverse momentum imbalance and the relative momentum of the quark pair. Obviously, the two variables are the conjugate momenta respect to v and u. In what follows, we are only interest in a special kinematic limit 2 , also known as the "back-to-back correlation limit" [7], in which the transverse momentum imbalance is much smaller than the relative momentum, i.e. |q ⊥ | |P ⊥ |. With this limitation, the corresponding variables u and u are much smaller than v and v , in this sense, terms involving multi-point Wilson line correlators can be expanded with Taylor series in coordinate space. To leading order in u and u 3 , Eq. (8) can be rewritten as, From the above expression, one can see that the differential cross section can be factorized into two parts, one is involved hard scale P ⊥ , which is referred to as hard part, the left one is related to q ⊥ , this is the quantity that directly probe the transverse momentum of incoming gluon from the hadronic target. Together with the saturation scale Q sg , there are three scales in this scattering process. Notice that, although |q ⊥ | |P ⊥ |, the value of q ⊥ can be order of (or softer than) the saturation scale Q sg , in this case, gluon saturation effects still remain important and this process can indeed provide information on the internal dynamics in the dense regime.
Within the back-to-back correlation limit, due to the two separated scales, i.e. |P ⊥ | |q ⊥ |, there is another class of effects with a double logarithmic form ln can give large contribution to the cross section and this should be included. This factor, also known as Sudakov resummations, comes from the resummation of soft collinear gluon radiation of incoming gluons from the scattered target. In the considered process, on the cross section level, which can be obtained by calculating one-loop correction of the quark pair [43,44]. After including this contribution, Eq. (9) becomes, To leading order, the Sudakov factor can be expressed as [44] where c 0 = 2e −γ E with γ E is the Euler constant.
It should be noted that the Sudakov factor only resums the soft radiations of incoming gluons from the scattered target, and this contribution is independent of u and u . Therefore, one can straightforwardly integrate over u and u by collecting Eq. (2) and Eq. (3), then Eq. (10) becomes, Due to P ⊥ is the conjugate momentum respect to u, the orientation of u i holds in the above expressions, and is represented by P i ⊥ after taking the Fourier transformation. Another piece is the double derivation of quadrupole, this one is related to the operator definition of WW gluon distribution. When averaging over impact parameter , it can be expressed as in which S ⊥ = d 2 B ⊥ is the transverse area of the scattered target. Generally speaking, xW ij is a 2 × 2 matrix, it can be decomposed into various different tensor structures. In what follows, it is convenient to separate it as trace and traceless parts, for positive definite, the corresponding tensor structure in coordinate space takes the form, where xG(x g , r ⊥ ) and xH(x g , r ⊥ ) are unpolarized and linearly polarized WW gluon distributions in coordinate space.
Combining Eq. (14) and Eq. (15), then the cross section becomes where ∆φ is the angle between P ⊥ and q ⊥ , dP.S. = d 3 k 1 d 3 k 2 , and xG(x g , P ⊥ , q ⊥ ) and xH(x g , P ⊥ , q ⊥ ) are the resumed WW gluon distributions in momentum space. From the above expressions, one can clearly see that the two distributions have resum small-x logarithms and Sudakov double logarithms simultaneously. We would like to point out that the resummations performed here is similar to that in Refs. [43,44], alternatively, a more general description of TMD evolutions has been done in [60], in which the two resummations can be performed without regarding in a specific scattering process in the CGC formalism. Notice that, for the linearly polarized WW gluon distribution, xH(x g , P ⊥ , q ⊥ ), one can not perform this transformation directly, this is mainly due to the tensor structure 1 2 δ ij − 2 In practice, when performing this transformation, a factor of − cos(2θ) will be generated, where θ is the angle between q ⊥ and r ⊥ . From Eq. (16) and Eq. (17), one can see that the differential cross section can be formulated in terms of xG and xH. Especially, the coefficients of cos(2∆φ) terms in the above expressions can provide us the direct information of linearly polarized WW gluon distributions, xH. If one integrate out φ, then the coefficients will disappear and the cross section only dependent on the unpolarized ones, xG. This indicates that the two gluon distributions can be extracted by measuring the total cross section and the coefficient of cos(2∆φ) respectively. As argued in Refs. [23,24], since we are working in coordinate space, ∆φ = φ u − φ v with φ u and φ v being the azimuthal angles of u and v, while P ⊥ and q ⊥ are the conjugate momenta of u and v, therefore, it is transparent to explore the properties of xH by measuring the asymmetry of P ⊥ and q ⊥ . Furthermore, attention should be paid is that, for transversely and longitudinally polarized photons, there exists a sign change for the coefficients of cos(2∆φ), as shown in Eq. (16) and Eq. (17). Especially, if Q ≈ 0 (for quasi-photon scattering process), the coefficient will become a negative one.

III. GENERAL PROPERTIES OF WW GLUON DISTRIBUTIONS IN MCLERRAN-VENUGOPALAN MODEL
In the last section, we have discussed the quark-antiquark pair (or dijet) cross section that can be formulated in terms of xG and xH, and pointed out that the cos(2∆φ) coefficients directly relate to the linearly polarized ones. In order to calculate the cross sections, the explicit expressions of xG and xH should be specified. In what follows, we do not include the Sudakov factor and will apply the original McLerran-Venugopalan (MV) model [54][55][56] to study the behaviors of the two gluon distributions. In this model, the color charge density ρ a (x + , x) conforms to a Gaussian distribution and satisfies the following relation, where µ 2 xg is the squared density of color charges per area. The above relation only hold when quantum evolution effects are not large, but x g is sufficiently small such that the high density regime can be accessed. In this sense, the color charge density ρ a (x + , x) is regarded as static source, while for small-x gluons, they are treated as classical gauge field A − a due to the density is very large and can be expressed in terms of g s ρ via the Yang-Mills equation, where G 0 (x) is the gluon propagator. Using Eq. (20), one can obtain the two-point function of gluon field A − a , with By introducing a dimensional function It is straightforward to show that the quark dipole reads, and the gluon dipole as, with where Q 2 s = 2C F α s g 2 s µ 2 xg and Q 2 sg = 2N c α s g 2 s µ 2 xg are the quark and gluon saturation scales, respectively. Λ QCD introduced here is the QCD infrared scale. It should be noted that the logarithmic term in the above expressions is very important and is the crucial difference between the Golec-Biernat-Wüsthoff (GBW) [61] and the MV model. Particularly, the logarithmic terms in MV model restores the correct high k ⊥ perturbative behaviors for the WW gluon distributions, ∼ Q 2 s /k 2 ⊥ (see for examples [62,63]). For the remainder of this work, we will use the functional form and choose Λ QCD = 0.241 GeV, Q 2 s0 = 0.2 GeV 2 as the inputs of rcBK equation. The two values used here are obtained by fitting the structure function F 2 with the HERA DIS ep data. We stress that the initial saturation scale Q 2 s0 is that of proton's at x 0 = 0.01, for heavy nucleus target by using Q 2 s0A = A 1/3 Q 2 s0 . Due to the fact that the saturation scale will increase with decreasing x g , in what follows, we use a model independent definition via the relation [64], In this sense, Q 2 s can be obtained by using the solution of rcBK equation, while, it is only the saturation scale in the fundamental representation, the gluon saturation scale can be obtained by using Q 2 sg = N c /C F Q 2 s . Another important ingredient is the double derivative of quadrupole correlator, which is directly relate to the definition of WW gluon distribution, Eq. (14), and can be expressed as [7,27] 1 with A. Unpolarized and linearly polarized WW gluon distributions in coordinate space In order to give a systematical study, we firstly express the two gluon distributions in coordinate space. By collecting Eq. (14) and performing the projection onto 1 2 δ ij and 1 2 (δ ij − 2 , one can directly obtain, and To investigate the general properties of unpolarized and linearly polarized WW gluon distributions in coordinate space, we typically distinguish two regions. For small dipole size limit, |r ⊥ | 1/Q sg , one can easily check that xG(x g , r ⊥ ) In the small r ⊥ limit, xG(x g , r ⊥ ) behaves as ln 1 |r ⊥ |ΛQCD and xH(x g , r ⊥ ) to be a constant at given x g , the perturbative behaviors are actually in agreement with the properties at large transverse momentum [11]. Moreover, it is interesting to note that the two gluon distributions have strong dependence on A since Q 2 sg ∼ A 1/3 and S ⊥ ∼ R 2 A ∼ A 2/3 , which indicates, in the small r ⊥ limit, the two gluon distributions have no nuclear shadowing effect, this feature also holds at large transverse momentum. Furthermore, it should be noted that |r ⊥ |ΛQCD , the ratio seems to be zero when r ⊥ is small enough even though xH still remains non-zero value. Last but not least, the evolutions of the two gluon distributions toward to small-x are completely captured by Q 2 sg . However, for , the two gluon distributions are independent on x g and fall off as r 2 ⊥ with increasing |r ⊥ |, these behaviors are completely different from that in small r ⊥ limit. Also, the ratio xH(xg,r ⊥ ) xG(xg,r ⊥ ) ≈ 1/2, which indicates the linearly polarized WW gluon distribution will tend to its saturated boundary 4 when r ⊥ is large enough.
In order to give numerical results of unpolarized and linearly polarized WW gluon distributions in coordinate space, according to Eq. (32) and Eq. (33), gluon saturation scale Q 2 sg should be known, this value can be obtained by solving the rcBK equation. To avoid repetition, details can be found in Refs. [48,50]. We stress that the key input of rcBK equation is the initial saturation scale Q 2 s0 , in what follows, we choose Q 2 s0 = 0.2 GeV 2 for proton. We firstly show the small-x evolutions of the two gluon distributions in coordinate space. In the left panel of Fig. 1, we present the numerical results of unpolarized and linearly polarized WW gluon distributions at different values of x g . Black curves denote the unpolarized WW gluon distributions and the red ones for the linearly polarized distributions. Solid, dashed and dot-dashed lines correspond to x g = 10 −2 , 10 −3 and 10 −4 , respectively. In the small r ⊥ region, from this figure, one can clearly see that the unpolarized WW gluon distributions warmly increase with decreasing r ⊥ , precisely behaves as ln 1 |r ⊥ |ΛQCD . While, for the linearly polarized ones, the results tend to flat over a broad small r ⊥ region, which indicates the linearly polarized WW gluon distributions still remains non-zero values even when r ⊥ is small enough. Besides that, one can observe that both of the two gluon distributions increasing with decreasing x g , the feature is expected since xG, xH ∼ Q 2 sg and the gluon saturation scale Q 2 sg goes large when x g becomes small. As in large r ⊥ region, e.g. |r ⊥ | ≥ 1 GeV −1 , one can observe that the unpolarized (linearly polarized) WW gluon distributions respective to different values of x g merge into a single curve, this is mainly due to the two gluon distributions are independent on x g when r ⊥ is large enough. Especially, in large r ⊥ region, one can see that both the two gluon distributions have same behaviors, precisely equal to r −2 ⊥ . In the right panel of Fig. 1, we present the ratio of xH(x g , r ⊥ )/xG(x g , r ⊥ ) as a function of r ⊥ . From this figure, one can observe that the numerical results for the three different values of x g merge into a single curve, this is 4 In practice, the saturated boundary is sensitive to the treatment of logarithmic term ln 1 |r ⊥ |Λ QCD . If one absorb 1 2 into the logarithmic term, then the saturated boundary tends to 1. mainly due to the ratio becomes a function of r ⊥ by comparing Eq. (32) and Eq. (33). Another feature is that the ratio goes large when increasing r ⊥ and tends to saturated boundary 1/2 when r ⊥ is large enough, which indicates xH(x g , r ⊥ ) ≤ 1 2 xG(x g , r ⊥ ). However, for physical aspect, xH(x g , r ⊥ ) ≤ xG(x g , r ⊥ ), this discrepancy mainly comes from the treatment of the logarithmic term in Eq. (33), ln −1 1 |r ⊥ |ΛQCD . In our calculations, we have ln −1 QCD + e , which will give the saturated boundary to be 1. We stress that the different treatments will do not alter the following analysis.

B. Unpolarized and linearly polarized WW gluon distributions in momentum space
Correspondingly, the unpolarized and linearly polarized WW gluon distributions can be obtained by using Fourier transformation, they are expressed as follows: and It is worthwhile to note that the linearly polarized WW gluon distribution cannot be transformed from coordinate space into momentum space directly, this is because the tensor structure 1 2 (δ ij − 2 ) is different from δ ij . After performing this projection, the tensor structure will naturally generate a − cos(2θ) factor, where θ is the relative angle between r ⊥ and k ⊥ . Alternatively, xG(x g , k ⊥ ) and xH(x g , k ⊥ ) can be arrived at by projecting xW ij onto 1 2 δ ij and , respectively, which has been done in Refs. [23,24].
In the above expression, J 2 (|k ⊥ ||r ⊥ |) is the Bessel function of the first kind. The following analysis is similar to that for gluon distributions in coordinate space. For k 2 ⊥ Q 2 sg , the dominated contribution comes from the small dipole size |r ⊥ | 1/Q sg , and can be evaluated by expanding the exponential term, then xG(x g , k ⊥ ) , in this limit, the two gluon distributions are completely same and behave as The perturbative property is agreement with other QCD results [62,63,65,66]. Moreover, due to xG( , the two gluon distributions scale like A and have no nuclear shadowing effects, therefore, each of nucleon can additively give a contribution to the distributions in the dilute region. This conclusion is also agreement with that in coordinate space at small r ⊥ limit. In addition, xH(xg,k ⊥ ) xG(xg,k ⊥ ) ≈ 1 when k ⊥ is large enough, which gives a saturated boundary for the linearly polarized WW gluon distribution in the dilute regime [23]. For Λ 2 QCD ≤ k 2 ⊥ Q 2 sg , the dominant contribution comes from large distances r ⊥ 1/Q sg , one can neglect the exponential and the logarithmic r ⊥ terms, and arrive at xG(x g , k ⊥ ) , one can check that the two gluon distributions scale like A 2/3 , which means the two gluon distributions have strong nuclear shadowing effects in the dense regime. Another interesting feature is the linearly polarized WW gluon distribution is weakly dependent on k ⊥ and approximates to a constant in this region. Furthermore, one can find xH(xg,k ⊥ ) xG(xg,k ⊥ ) is suppressed in the small k ⊥ limit, in which the multi-gluon scattering effects play a role that constrains the gluons to be polarized. In Fig. 2, we present the numerical results of the two gluon distributions in momentum space. In the left panel of Fig. 2, black and red curves denote the unpolarized and linearly polarized WW gluon distributions, solid, dashed and dot-dashed lines correspond to the results at x g = 10 −2 , 10 −3 and 10 −4 , respectively. As shown in the left panel of Fig. 2, for each value of x g , one can clearly see that the results for unpolarized and linearly polarized WW gluon distributions approach to a single line when k ⊥ is very large, which implies the linearly polarized WW gluon distribution tends to its saturated boundary, xH(x g , k ⊥ ) ≈ xG(x g , k ⊥ ). Another important feature is that both the two distributions have same asymptotic behavior, i.e. an inverse power law, precisely equal to k −2 ⊥ when k ⊥ lies in large transverse momentum region. Furthermore, with decreasing x g , one can observe the curves are shifted to large k ⊥ , this is because, in the large k ⊥ region, the evolutions of the two gluon distributions toward to small x are fully characterized by gluon saturation scale, Q 2 sg , and which is the value becomes large when x g decreases. As k ⊥ goes small, the two gluon distributions start to separate. Specifically, the unpolarized gluon distributions will warmly increase with decreasing k ⊥ , while the linearly polarized ones tend to flat over a broad range of k ⊥ and are almost independent on x g , this feature depicted here is again confirms what we discussed in the above.
In the right panel of Fig. 2, we present the ratio of xH(x g , k ⊥ )/xG(x g , k ⊥ ) as a function of |k ⊥ |/Q sg for three different values of x g . The solid, dashed and dot-dashed lines correspond to the results at x g = 10 −2 , 10 −3 and 10 −4 , respectively. From this figure, one can see that these quantities are weakly dependent on x g and approx to 1 when |k ⊥ |/Q sg > 4 ÷ 5, this strongly indicates the geometric scaling behavior [67,68] is held for the two gluon distributions and the saturated boundary of linearly polarized WW gluon distribution is reached at large k ⊥ region. When k ⊥ goes small, such as of order of Q sg , in which gluon multiple rescattering (saturation) effects to be manifest and play a role that suppressing the gluons to be polarized. Therefore, as shown in this figure, the effects will naturally lead the ratio to be smaller. We note that, although the ratio goes small when decreasing k ⊥ , it is still sizable even for small k ⊥ values, xH(x g , k ⊥ )/xG(x g , k ⊥ ) ≥ 10%. In this sense, the dijet cos(2∆φ) structure 5 in the considered process is expected to be visible even though in small transverse momentum region. For completeness, we also present small-x evolutions for the two gluon distributions, results are shown in Fig. 3. In the calculations, we have chosen two values of k ⊥ , results presented in the left panel of Fig. 3 correspond to the case of |k ⊥ | = 1 GeV, and the right panel is for |k ⊥ | = 2 GeV. The solid (black) and dashed (red) curves denote unpolarized and linearly polarized WW gluon distributions, and the dot-dashed (blue) vertical lines represent the values of x g by setting Q sg = |k ⊥ |, respectively. In the left panel of Fig. 3, one can see that the vertical line separates the domain as two regions, specifically, in which x g ≥ 0.003 with |k ⊥ | ≥ Q sg as the dilute regime and x g ≤ 0.003 with |k ⊥ | ≤ Q sg as the dense (saturated) regime. As discussed in the above, we recall that the gluon saturation effects will be very significant in the dense regime (x g ≤ 0.003), where gluon rescattering effects will be very important and constrain those gluons to be polarized, therefore, the linearly polarized part becomes flat as shown in the figure. While, in the dilute regime (x g ≥ 0.003), the effect is less manifest and makes the linearly polarized part close to the unpolarized one. This feature is more visible for the case of |k ⊥ | = 2 GeV (see the right panel of Fig. 3).

IV. NUMERICAL RESULTS AND DISCUSSIONS
According to Eq. (16) and Eq. (17), the quark-antiquark pair (or dijet) cross section is directly related to unpolarized and linearly polarized WW gluon distributions, when the two distributions are known, then the cross section can be obtained. To investigate the azimuthal correlations of dijet in γ * p(A) scattering processes, we define the azimuthal correlation function in the following form, where ∆φ is the difference between the azimuthal angles of P ⊥ and q ⊥ . This definition has a clear physical interpretation that quantifying the probability distribution of the dijet versus the relative angle ∆φ.
To perform numerical calculations, we will focus on typical EIC kinematics and take the center-of-mass energy of the scattering system to be W = 140 GeV. When the kinematics are specified, x g enters Eq. (16) and Eq. (17) and can be determined via, where M p is the mass of scattering proton. When performing calculations, we also restrict |P ⊥ | ≥ 2|q ⊥ | in order to meet the back-to-back correlation limit. We stress that the chosen kinematics in the following is also relevant for probing gluon saturation effects of scattered target, as we will demonstrate below. In addition, according to Eq. (37), x g increases with increasing P ⊥ , while it has weak dependence on q ⊥ , which indicates that dσ/d|P ⊥ | can be used as a proxy to measure the evolutions of xG(x g , q ⊥ ) and xH(x g , q ⊥ ) towards to x g , while dσ/d|q ⊥ | as a direct probe that measuring the transverse momentum of gluon inside the scattered target. In Fig. 4, we show the differential cross section for longitudinally and transversely polarized photons as a function of transverse momentum |P ⊥ |, in which we have set the transverse momentum imbalance |q ⊥ | = 0.5 GeV and restrict the dijet to be purely back-to-back, ∆φ = π. For the phenomenological study, we also vary the virtuality of photons as order of transverse momentum, specifically, the left panel in Fig. 4 shows the results for the case of Q = 2|P ⊥ | and the right panel corresponds to the case of Q = |P ⊥ |. It also should be reminded that the presented results are for light flavor quark jets with quark mass m q = 0.01 GeV and symmetric longitudinal jet z = 1/2. As can be seen from this figure, the cross sections are rapidly falling with increasing the dijet transverse momentum |P ⊥ |. For different virtualities Q, the contribution of longitudinally photons and transversely photons to the total cross section is different. When Q = 2|P ⊥ | (see left panel of Fig. 4), one can observe that the longitudinal part is larger than the transverse part, in this case, the total cross section is mainly driven by the longitudinal part. For another case, when Q = |P ⊥ | (see right panel of Fig. 4), although the longitudinal part is still larger than the transverse part, the two parts seem to be comparable in the total cross section. By comparing the two figures, one can conclude that the contribution of transverse part to the total cross section will increase with decreasing the photon's virtualities, Q 2 . Similar to Fig. 4, in Fig. 5 we show the differential cross section as a function of transverse momentum imbalance |q ⊥ |, in which we have set the transverse momentum |P ⊥ | = 4 GeV and restrict the dijet to be purely back-to-back, ∆φ = π, the other variables are same as in Fig. 4. From this figure, one can observe that the contributions of transversely and longitudinally polarized photons to the total cross section are quite different as comparing to that of in Fig. 4. Specifically, for the case of Q = 2|P ⊥ | (left panel of Fig. 5), one can observe the distribution of transverse part decreases much more sharply than the longitudinal part as increasing q ⊥ , therefore, the total cross section is almost driven by the longitudinally polarized photons when q ⊥ is enough large. The peculiar behaviors exhibited by transversely and longitudinally photons in the figure are not surprising actually, and can be understood as follows: Before to clarify that, it is useful to remind that the transverse momentum imbalance q ⊥ is the only variable in the calculations, and as the quantity that directly probes the transverse momentum of incoming gluons from the scattered target. According to Eq. (16) and Eq. (17), if only the unpolarized WW gluon distribution xG(x g , q ⊥ ) is considered or integrating over ∆φ, one can expect that the two cross sections are parallelly falling with increasing q ⊥ . Moreover, it should be noted that the dijets are restricted to be purely back-to-back (∆φ = π), in this sense, the cos(2∆φ) terms will give negative (positive) contribution respect for transversely Eq. (16) (longitudinally Eq. (17)) polarized photons, therefore, the two cross sections will naturally be separated and the separation between the two is directly related to the quantity of xH(x g , q ⊥ ) or the ratio of xH(x g , q ⊥ )/xG(x g , q ⊥ ). Furthermore, in the specific kinematics, since x g is weakly dependent on q ⊥ and x g ≈ 0.0065, correspondingly, the gluon saturation scale Q sg ≈ 0.88 GeV as shown by the vertical line in the left panel of Fig. 5. Reminding that the ratio of xH(x g , q ⊥ )/xG(x g , q ⊥ ) becomes large when q ⊥ goes large (also see the right panel of Fig. 2). As a consequence, the separation between the two cross sections will become more significant as |q ⊥ | ≥ Q sg , which is shown in the left panel of Fig. 5. Similar trends hold for the case of Q = |P ⊥ | (right panel of Fig. 5) as well.
As mentioned in previous, the longitudinally polarized photons (Eq. (17)) will give a cos(2∆φ) structure in azimuthal correlations, in contrast, the transversely polarized photons (Eq. (16)) gives a − cos(2∆φ) structure 6 , and this part plays a role that partially washing out the cos(2∆φ) structure in azimuthal correlations. One can expect that the cos(2∆φ) structure will become more visible only when the total cross section is dominated by longitudinally polarized photons, otherwise, if the total cross section is mainly driven by the transverse part, it is a challenge to observe the cos(2∆φ) structure. In Fig. 6, we plot the normalized azimuthal correlations C(∆φ) defined in Eq. (36) as a function of relative angle ∆φ. In the left panel of Fig. 6, the represented results are for different photon virtualities Q by setting |P ⊥ | = 2 GeV and |q ⊥ | = 0.5 GeV, the other kinematics are same as in Fig. 4. The solid (black), dashed (red) and dot-dashed (blue) lines correspond to the case of Q = 2|P ⊥ |, |P ⊥ | and 1 GeV, respectively. As can be seen from this figure, the cos(2∆φ) structure in azimuthal correlations will be suppressed and broaden with decreasing the virtualities Q, this is mainly due to the fact that the contribution of longitudinally polarized photons to the total cross section will decrease as Q decreases. In the right panel of Fig. 6, we show the azimuthal correlations for different q ⊥ , in the calculations, we set Q = 2|P ⊥ | with |P ⊥ | = 2 GeV, the other kinematics are same as former. The solid (black) and dashed (red) lines correspond to |q ⊥ | = 0.5 and 1 GeV, respectively. As discussed in previous, the cos(2∆φ) structure in azimuthal correlation is sensitive to the ratio of xH(x g , q ⊥ )/xG(x g , q ⊥ ), and the ratio will become large as q ⊥ goes large, therefore, one can observe that the cos(2∆φ) structure will become more visible for |q ⊥ | = 1 GeV as compared to that of |q ⊥ | = 0.5 GeV.
Furthermore, we have also checked that for the heavy flavor jets, where the heavy quark mass takes m q = 1.4 GeV, results are very similar to the light quark jets. For this reason, we do not include separate plots.

B. Sudakov factor will suppress the azimuthal correlations
In order to further investigate how does the Sudakov resummations impact on the cross sections, we perform calculations by regarding the same kinematics as in Fig. 4, results are shown in Fig. 7. Similarly, the left and right panels in Fig. 7 correspond to the case of Q = 2|P ⊥ | and Q = |P ⊥ |, the black and red curves denote without and with Sudakov resummations, respectively. By comparing the black and red curves depicted in Fig. 7, one can observe that the behavior of the cross sections is not changed after including the Sudakov resummations, the cross sections are tamely falling with increasing |P ⊥ |. However, a difference between the two (black and red curves) is that the Sudakov resummations will suppress the cross sections and the suppression tends to be significant when |P ⊥ | goes large, as revealed by the magnitude in the figure. It naturally reflects the fact that the parton shower effect is relatively weaker in the low |P ⊥ | region while the effect goes stronger in the large |P ⊥ | region. This argument is consistent with other perturbative calculations, see e.g. Ref. [41]. Another interesting feature by comparing the black and red curves, such as for the case of Q = 2|P ⊥ | (left panel of Fig. 7), is that the relative difference between transverse and longitudinal parts goes smaller when the Sudakov resummations are taken into account. More precisely, for the specific kinematics, the ratio of transverse component in the total cross section σ T /σ Total ≈ 0.32 when including the Sudakov resummations, while σ T /σ Total ≈ 0.25 if the Sudakov resummations are not included. As discussed in previous, transversely polarized photons gives a contribution that partially washing out the cos(2∆φ) structure in azimuthal correlations, therefore, one can straightforwardly expect that the cos(2∆φ) structure will be suppressed and broaden if the Sudakov resummations are taken into account. We further notice that, after including Sudakov resummations, the transverse part seems to be parallelly falling with increasing P ⊥ respect to the longitudinal part, in this sense, the cos(2∆φ) structure in azimuthal correlations will not be changed for different values of P ⊥ . Similar features hold for the case of Q = |P ⊥ | (see right panel of Fig. 7). In Fig. 8, we show the cross sections as a function of transverse momentum imbalance |q ⊥ |, the other kinematics are same as in Fig. 5. Similarly, the black and red curves denote without and with Sudakov resummations, the left and right panels correspond to the case of Q = 2|P ⊥ | and Q = |P ⊥ |, respectively. As shown in this figure, the red curves are significantly suppressed as compared to the black one in small q ⊥ region, while with increasing q ⊥ , one can observe this suppression goes small 7 . This actually reflects the fact that Sudakov effects will reduce the distribution at small q ⊥ region and make the distributions to be flattened over a broad q ⊥ region. We notice that the distributions presented here are somewhat consistent with a recent study in Ref. [69], in which the authors have pointed out the Sudakov resummations are very important in describing the dijets transverse momentum imbalance distributions at E791 [70,71], and shown only the Sudakov effects is considered that can make the theoretical results meet to the experimental data. Another important feature as shown in the figure is that the difference between transverse and longitudinal cross sections will become smaller if the Sudakov resummations are considered, therefore, the cos(2∆φ) structure in azimuthal correlations is expected to be suppressed and broaden. In Fig. 9, we show the normalized azimuthal correlations C(∆φ) as a function of ∆φ. In the left panel of Fig. 9, we present the results for different P ⊥ , in which we have set Q = 2|P ⊥ |, |q ⊥ | = 0.5 GeV and z = 1/2. The solid, dashed and dot-dashed lines correspond to |P ⊥ | = 4, 2 and 1 GeV, the red and black curves denote with and without Sudakov resummations, respectively. From this figure, one can observe that the back-to-back (∆φ = π) correlations are indeed suppressed when including Sudakov resummations, especially, this suppression goes stronger when P ⊥ goes larger. Moreover, in the specific kinematics, one can observe that the red curves are completely overlap together, which indicates the cos(2∆φ) structure in azimuthal correlations will be unchanged after including Sudakov resummations. This feature can be tracked back to that the transverse cross sections will be parallelly falling respect to the longitudinal ones as P ⊥ goes large (see for example Fig. 7).
As mentioned in previous, the effect of Sudakov resummations is strongly sensitive to the difference between |P ⊥ | and |q ⊥ |, when Q and |P ⊥ | are fixed, one may expect that the effect will decrease with increasing |q ⊥ |. In order to verify this, we present results in the right panel of Fig. 9, where Q = 4 GeV and |P ⊥ | = 2 GeV, the solid and dashed curves correspond to |q ⊥ | = 0.5 and 1 GeV, the black and red curves denote without and with Sudakov resummations, respectively. From this figure, besides the fact that the Sudakov resummations suppressing the backto-back correlations, one can observe that the suppression for |q ⊥ | = 0.5 GeV is stronger than that of |q ⊥ | = 1 GeV (By comparing the differences between the away-side peaks for the two cases.). Furthermore, it is worthwhile to note that, after including Sudakov resummations, the away-side (∆φ = π) peak for |q ⊥ | = 1 GeV is still larger than that of |q ⊥ | = 0.5 GeV, the main reason comes from the fact that the cos(2∆φ) structure is sensitive to xH(x g , q ⊥ )/xG(x g , q ⊥ ) and the ratio becomes large as q ⊥ goes large.
C. Gluon saturation further suppresses the azimuthal correlations of the dijets As mentioned in previous, the cos(2∆φ) terms in cross sections directly quantify the values of linearly polarized WW gluon distributions (Eq. (16) and Eq. (17)), this distribution is strongly suppressed when the gluon saturation dominated.
To further demonstrate the gluon saturation effects will actually suppress the back-to-back correlations, we plan to perform a calculation to compare the dijets azimuthal correlations in γ * p and γ * A scattering systems with the same kinematics, in which the only difference is the initial saturation scale, for γ * p scattering process, the initial saturation scale for proton Q 2 s0p = 0.2 GeV 2 , for γ * A, we choose Q 2 s0A = 1.2 GeV 2 due to Q 2 sA ∼ A 1/3 Q 2 sp with A ≈ 200. Within this assumption, one can expect that the back-to-back correlation in γ * A scattering process will be suppressed as compared to that in γ * p if gluon saturation effects is the actual reason that smearing the correlations. In Fig. 10, we present the normalized azimuthal correlations for different scattering systems. The black and red curves denote γ * p and γ * A scattering processes, the solid and dashed lines correspond to without and with Sudakov resummations. The left panel is for Q = 2|P ⊥ | and the right panel is for Q = |P ⊥ |, in which we have set |P ⊥ | = 2 GeV, |q ⊥ | = 0.5 and z = 1/2. From this figure, by comparing the solid lines with dashed lines, one observe that the back-to-back (∆φ = π) correlations are suppressed for each colliding system, as discussed in previous, this suppression is mainly come from Sudakov resummations. It should be noted that the chosen kinematics for γ * p and γ * A systems are same in the calculations, therefore, the effects of Sudakov resummations should be same in the two scattering systems. Furthermore, by comparing the red curves with black ones, one can see that the dijets back-toback (∆φ = π) correlations in γ * A are significantly suppressed than that of γ * p. This suppression is mainly due to the gluon saturation effects for heavy ion A is more predominant than the proton p. To demonstrate that, it is useful to remind that the enhancement of back-to-back correlations is proportional to the ratio of |q ⊥ |/Q sg , and with the fact of gluon saturation scale of A is larger than that of proton (Q 2 sA = A 1/3 Q 2 sp ), therefore, when the kinematics are fixed, the values of |q ⊥ |/Q sg in γ * A scattering system is smaller than that in γ * p, as a consequence, the back-to-back correlations in γ * A is naturally suppressed as compared to that in γ * p scattering process.
In addition, we also investigate the cos(2∆φ) anisotropy in the two scattering systems, which is the quantity that can be directly measured in experiments and is defined as follows, v 2 ≡ d∆φ cos(2∆φ)dσ d∆φdσ (38) Clearly, the numerator d∆φ cos(2∆φ)dσ ∼ xH(x g , q ⊥ ), while the denominator d∆φdσ ∼ xG(x g , q ⊥ ), hence, v 2 ∼ xH(x g , q ⊥ )/xG(x g , q ⊥ ). Therefore, by measuring the elliptic anisotropy v 2 of the dijets, we can get the information about xH/xG in the scattered target. In Fig. 11, we present the results of elliptic anisotropy v 2 as a function of transverse momentum imbalance q ⊥ for different scattering systems, in which we have set |P ⊥ | = 4 GeV and z = 1/2. The black and red curves correspond to γ * p and γ * A processes, the solid and dashed lines denote without and with Sudakov resummations, respectively. The left panel corresponds to the case of Q = 2|P ⊥ | and the right panel is for Q = |P ⊥ |. From this figure, for example in the left panel of Fig. 11, one can observe that the elliptic anisotropy v 2 for γ * A is significantly lower than that of γ * p, the main reason comes from the fact that the gluon saturation effects for heavy ion A is stronger than that of proton p. By comparing the solid lines with dashed ones, one can see that the soft gluon resummations significantly reduces the elliptic anisotropy v 2 as well. Especially, in the present kinematic regime, we stress that, for γ * A, the rather strong suppression of v 2 caused by the effects of gluon saturation and soft gluon resummations will make it very difficult to measure in experiments. On the other hand, the strong suppression for γ * A as compared to that of γ * p also as a signature of gluon saturation for heavy nucleus ion. In order to give a sizable elliptic anisotropy v 2 that can be measured in experiments, one possible way is to enlarge q ⊥ or increase x g that making |q ⊥ |/Q sg not even small.
Although the elliptic anisotropy v 2 can be used to probe xH/xG of the scattered target, it involves the other hard factors, such as Q, P ⊥ , etc. Following the same manner, the polarized anisotropy can be expressed [21,22], and v T 2 ≡ d∆φ cos(2∆φ)dσ T d∆φdσ T = − 1 2 R xH xG (40) where R is determined by kinematics in Eq. (16) From Eq. (39) and Eq. (40), one can observe that the longitudinally polarized elliptic anisotropy v L 2 directly probes the ratio of xH/xG and has no dependence on other additional variables, while for the transversely polarized anisotropy v T 2 , which involves a factor of R and has a sign change as compared to v L 2 . One should note that these polarized anisotropies can not be directly measured in experiments, one possible way to extract the polarized anisotropies following the relation where R is determined by comparing Eq. (16) and Eq. (17) In addition, we have also computed the longitudinally polarized anisotropy v L 2 in γ * p and γ * A scattering processes, the results are very similar to that of v 2 , therefore, we do not include it as a separated plot.

V. SUMMARY
In this paper, we have investigated the inclusive dijet cross sections and their azimuthal correlations in ep and eA collisions within the framework of CGC EFT. By applying the so-called "back-to-back correlation limit", the cross section can be expressed in terms of unpolarized and linearly polarized WW gluon distributions. Although the unpolarized and linearly polarized WW gluon distributions have little known in experiments, the novel scattering processes provide the cleanest environment to probe the two gluon TMDs, as well as allow one to quantitatively probe the gluon saturation effect inside the scattered target.
As compared to previous related studies [21,22,27,47], in this work, we employed the non-linear Gaussian approximation, using the solutions of rcBK equation to give the numerical results of unpolarized and linearly polarized WW gluon distributions both in coordinate and momentum space. Results shown that the two gluon distributions have strong geometrical scaling behaviors over a broad range of transverse momentum. Especially, for the linearly polarized WW gluon distribution, it tends to be saturated and has same asymptotic behavior as the unpolarized one at large transverse momentum region, while in small transverse momentum region, in which gluon saturation effects become manifest, this distribution remains non-zero values even at small transverse momentum, and seems to be constant when transverse momentum is order of gluon saturation scale (see for example Fig. 2 and Fig. 3). We stress that, although only the two gluon TMDs are presented in the work, the other gluon TMDs such as enter the cross sections at pp and pA collisions [27,47,72] can also be obtained following this procedure, and hope these features depicted in the present work can be further validated by experiments in future.
Moreover, we performed the calculations of quark-antiquark pair (or dijet) cross sections and their azimuthal correlations by including the parton shower effects (or Sudakov resummations) in the scattering processes. In the calculations, we have restricted the transverse momentum P ⊥ and the transverse momentum imbalance q ⊥ of the dijet with a condition of |P ⊥ | ≥ 2|q ⊥ | to meet the "back-to-back correlation limit". For the cross section, although the dijet are restricted to be purely back-to-back (∆φ = π), we found that the Sudakov resummations can strongly affect the transverse momentum distributions. Specifically, for dσ/d|P ⊥ |, numerical results shown that the Sudakov resummations suppresses the cross sections and the suppression goes stronger with increasing P ⊥ (shown in Fig. 7), while for dσ/d|q ⊥ |, the Sudakov effects will significantly reduce the distributions as q ⊥ goes small, and seems to make the distributions to be flattened over a broad small q ⊥ region (Fig. 8). As for the dijet azimuthal correlations, numerical results shown that the back-to-back correlations are significantly suppressed when the Sudakov effects are taken into account (Fig. 9).
Furthermore, we also presented the results of dijet azimuthal correlations and elliptic anisotropy v 2 in the γ * p and γ * A scattering processes. By comparing the back-to-back correlations in the two scattering processes, we found that a significant suppression that caused by gluon saturation effects is visible, especially, when including the Sudakov effects, that make the back-to-back correlations to be more suppressed. We also found that similar suppression can be observed in the elliptic anisotropy v 2 as well. Especially, in the studied kinematic regime, we stress that, for γ * A, the rather strong suppression of v 2 caused by the effects of gluon saturation and soft gluon resummations will make it very difficult to measure for experiments as compared to that of γ * p (see Fig. 11). On the other hand, as compared to γ * p, the rather strong suppressions of back-to-back correlations and elliptic anisotropy in γ * A is a solid signature of gluon saturation for heavy nucleus, which can be further validated or tested at a future EIC.
We should note that, for the dijet production, although our numerical results are performed only on parton level, similar trends are expected to hold for dihadrons, this can be done by convoluting the dijet cross section with fragmentation functions. We leave it as future study.