Comparison of improved TMD and CGC frameworks in forward quark dijet production

For studying small-$x$ gluon saturation in forward dijet production in high-energy dilute-dense collisions, the improved TMD (ITMD) factorization formula was recently proposed. In the Color Glass Condensate (CGC) framework, it represents the leading term of an expansion in inverse powers of the hard scale. It contains the leading-twist TMD factorization formula relevant for small gluon's transverse momentum $k_t$, but also incorporates an all-order resummation of kinematical twists, resulting in a proper matching to high-energy factorization at large $k_t$. In this paper, we evaluate the accuracy of the ITMD formula quantitatively, for the case of quark dijet production in high-energy proton-proton($p+p$) and proton-nucleus ($p+A$) collisions at LHC energies. We do so by comparing the quark-antiquark azimuthal angle $\Delta\phi$ distribution to that obtained with the CGC formula. For a dijet with each quark momentum $p_t$ much larger than the target saturation scale, $Q_s$, the ITMD formula is a good approximation to the CGC formula in a wide range of azimuthal angle. It becomes less accurate as the jet $p_t$'s are lowered, as expected, due to the presence of genuine higher-twists contributions in the CGC framework, which represent multi-body scattering effects absent in the ITMD formula. We find that, as the hard jet momenta are lowered, the accuracy of ITMD start by deteriorating at small angles, in the high-energy-factorization regime, while in the TMD regime near $\Delta\phi=\pi$, very low values of $p_t$ are needed to see differences between the CGC and the ITMD formula. In addition, the genuine twists corrections to ITMD become visible for higher values of $p_t$ in $p+A$ collisions, compared to $p+p$ collisions, signaling that they are enhanced by the target saturation scale.


Introduction
Parton saturation at small Bjorken's x in hadron wave functions is one of the most salient and universal features of QCD dynamics [1][2][3]. Small-x partons are interpreted as short-lived quantum fluctuations splitting from larger-x partons in a hadron wave function. Lorentz time-dilation dictates that the higher the collision energy is, the smaller-x partons come to participate in the interaction. The x-evolution of the gluon density has been formulated as the Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) equation [4][5][6][7][8][9], or the Balitsky-Kovchegov (BK) equation [10,11] in a mean-field approximation. The evolution changes from a linear to a non-linear character when the gluon density becomes so dense that the gluon merging starts to compete with the splitting. This transition is characterized by the so-called saturation momentum scale, Q s (x) [1][2][3], an emergent scale in QCD dynamics. Then, the color-glass-condensate (CGC) framework [12][13][14][15], which describes the small-x part of the wave function in the presence of large-x random color sources, has been realized as a suitable effective theory to calculate observables in the dense gluon regime with Q s (x) Λ QCD . Forward dijet production in proton-proton (p + p) and proton-nucleus (p + A) collisions at the large hadron collider (LHC) is a unique and valuable observable among others for the phenomenological study of gluon saturation. In this process a large-x parton from the projectile, which is dilute and well understood in perturbative QCD, probes the small-x partons in the dense target, and then produces jets at forward rapidities. This setup is sometimes called dilute-dense system. In addition to its ever highest collision energy, the nuclear target option available at the LHC is very advantageous since gluon saturation, or its scale Q s (x), is enhanced by the target thickness ∝ A 1/3 (A is the nuclear mass number).
In the CGC framework, the dijet production cross-section is expressed in terms of the Wilson line correlators averaged over external color source distributions. The Wilson-line correlators with fixed transverse positions are essential components to define the gaugeinvariant matrix elements. Those correlators encode multiple scatterings of the partons traversing the dense target and satisfy the BK-JIMWLK evolution, provided that leading logarithms in x are predominant over leading logarithms in Q 2 . Those multiple scattering effects are enhanced in the dense regime where the saturation scale Q s increases. It is demonstrated in Ref. [16,17] that the description of dijet production at the LHC should simplify thanks to the hard scales involved there.
Indeed, the dijet production contains three characteristic momentum scales: the typical transverse momentum of a hard jet P t , the transverse momentum imbalance of the pair k t , and the saturation scale of the target Q s . Here P t is always the hardest scale, while Q s is the softest of the three. The original CGC framework does not assume any ordering in the three momentum scales. In the Q s P t ∼ k t limit, expanding the Wilson line correlators in the CGC expression to the second order in the gluon field, one can obtain the "dilute" result known as high-energy factorization (HEF) or k t -factorization. On the other hand, in the Q s ∼ k t P t limit, by keeping the leading 1/P t terms from the CGC expression, one can accurately reproduce the leading-twist TMD factorization result at small x which comes with on-shell hard matrix elements.
In the meantime, by introducing the off-shell k t dependence of the small-x gluons in the hard matrix elements, Ref. [16,17] proposed an improved TMD (ITMD) expression, which is valid for any k t provided Q s P t , and interpolates the TMD and HEF expressions. Then it was pointed out in Refs. [18,19] that such off-shell effect results from the resummation of power corrections in |k t |/|P t | in the hard scattering parts, known as kinematic-twists corrections, coupled to leading-twist TMD distributions. Alternatively, the ITMD framework can also be thought of as an improvement of HEF, from that perspective the HEF framework gets supplemented with leading-twist saturation corrections. The ITMD framework provides a concise and useful approximation to the CGC expression for Q s P t , and it is crucial now to assess the quantitative accuracy of the ITMD formula, compared to the "full" CGC formula, when calculating the spectrum of forward dijets. This is a practical motivation of this paper.
Gluon saturation affects particle production in hadron collisions through the non-linear evolution of the gluon density, and through the multiple scattering of the partons with the dense target. The multiple scattering effects are further categorized into two classes: the leading-twist ones accounted for in the (I)TMD framework, controlled by the magnitude of |k t | vs. Q s , and those due to genuine higher-twist effects, controlled by |P t | vs. Q s . The CGC formula contains both effects of multiple scatterings, while the ITMD formula is obtained from the CGC one by getting rid of the genuine higher-twist corrections, which may be referred to as Wandzura-Wilczeck approximation [20]. The numerical comparison of the ITMD to the CGC formula will give valuable information about the genuine highertwist effects on forward dijets production in high-energy p + A collisions. In order to make our ITMD/CGC comparison feasible and clear, we shall restrict our analysis to the forward quark (qq) dijet production, and work within the Gaussian truncation of JIMWLK evolution and large-N c limit, for which the CGC expression is less complicated and can be evaluated directly (indeed, as we will see below, the two expressions then differ only in their hard factors). In this regard, we note that genuine-twist corrections were also analyzed recently in the context of dijet production in deep-inelastic scattering [21], using the same approximation but keeping finite N c corrections.
The paper is organized as follows; Section 2 gives an overview of the ITMD and CGC frameworks for forward dijet production. In Section 3, we present numerical results on the dijet azimuthal angle correlation in the ITMD and CGC frameworks. In particular, we will look into the dependence of the genuine-twist corrections on kinematics and system size there. Section 4 is devoted to summary and concluding remarks.

Frameworks
This section runs through some details of the ITMD and CGC frameworks for forward dijet production in dilute-dense collisions.

Improved TMD factorization for forward dijet production
We consider the process of inclusive dijet production at forward rapidity in proton-nucleus collisions where the four-momenta of the projectile and the target are massless and purely longitudinal. In terms of the light cone variables, x ± = (x 0 ± x 3 )/ √ 2, they take the simple form p µ p = s/2 (1, 0, 0 t ) and p µ A = s/2 (0, 1, 0 t ) where s is the squared center of mass energy (per nucleon-nucleon collisions) of the p + A system. The longitudinal momentum fractions of the incoming parton from the projectile, x 1 , and of the gluon from the target, x 2 , can be expressed in terms of the rapidities (y 1 , y 2 ) and transverse momenta (p 1t , p 2t ) of the produced jets as By looking at jets produced in the forward direction, we effectively select those fractions to be x 1 ∼ 1 and x 2 1. Since the target A is probed at low x 2 , the dominant contributions come from the subprocesses in which the incoming parton on the target side is a gluon, meaning there are three possible channels: qg → qg, gg → qq, and gg → gg. Figure 1 shows the kinematics for the gg → qq subprocess in p + A collisions.
The asymmetry of the problem, x 1 ∼ 1 and x 2 1, also implies that gluons from the target have a much bigger average transverse momentum (of the order of Q s (x 2 )) compared to that of the partons from the projectile (which is of the order of Λ QCD ). Therefore we shall always neglect the transverse momentum of the high-x 1 partons from the projectile compared to that of the low-x 2 gluons from the target. As a result, the parton content of the projectile hadron is described by regular collinear parton distributions f a/p (x 1 , µ 2 ) (where µ is the factorization scale) and TMDs are involved only on the target side, with the transverse momentum of those small-x 2 gluons being equal to the transverse momentum of jet pair k t : This simplification is needed to apply the TMD factorization for the dijet process, since for this final state, there is no such factorization with TMDs for both incoming hadrons [22,23]. The ITMD factorization formula reads [16] dσ(p where P t is the hard scale of the process, related to the individual jet momenta: with z = p + 1 /(p + 1 + p + 2 ) the longitudinal momentum fraction carried by the jet j 1 . The improvement with respect to the TMD factorization formula derived in Ref. [24] (in the large-N c limit) and in Ref. [25] (keeping N c finite), lies in the fact that the hard factors H (i) ag * →cd (P t , k t ) are k t -dependent, as opposed to a function of P t only in the TMD case: H ag * →cd (P t , 0); their expressions can be found in Ref. [16]. On the other hand, the improvement with respect to the HEF lies in the fact that several gluon distributions are involved, which differ from one another when non-linear effects become important. The various operator definitions of the gluon TMDs F (i) ag (x 2 , k t ) are found in Ref. [25]. From now on, we focus solely on a quark dijet pair (qq) production, since considering this subprocess will allow us to make a detailed comparison with the CGC formulation. In that case, let us write down more explicitly the ITMD formula 1 ; 1 Compared to Ref. [25], F gg is simply denoted Fgg, the Weizsäcker-Williams gluon TMD F gg is denoted FW W , and F adj = F (1) gg − F (2) gg is the adjoint-dipole gluon TMD [26].
x 1 where denotes the usual gluon-quark splitting function at leading order in α s . The relevant small-x gluon TMDs are given by [25] F gg (x 2 , k t ) = 4 g 2 in terms of the Wilson lines with t a and T a denoting the generators of the fundamental and adjoint representation of SU (N c ), respectively. The process is depicted in Fig. 1 at the amplitude level, the ITMD cross-section being the square of Fig. 1 (a). The soft gluons attaching to the hard parts, are not shown, those are accounted for by two (fundamental) Wilson lines. A derivative applied to a Wilson line corresponds to a gluon exchanged in the t-channel, those are explicitly drawn. The CGC averages · x 2 represent averages over the configurations of the classical color field of the hadronic/nuclear target, A − , which describes the dense parton content of its wave function, at small longitudinal momentum fraction x 2 . In the leading-logarithmic approximation, the evolution of the CGC averages with decreasing x 2 obeys the JIMWLK equation, where H JIMWLK denotes the JIMWLK Hamiltonian.
The ITMD formula (2.4) is an interpolation between two limiting cases, Q s |k t |, |P t | and Q s , |k t | |P t |, both limits being contained as well in the more general CGC framework (the details of which are recalled below). The ITMD formula is valid when |P t | Q s (x 2 ), however the value of k t can be arbitrary. When k t Q s (x 2 ), the high-energy factorization (HEF) formula (aka k t -factorization) is recovered: the various gluon TMDs (2.8) collapse into a single function, known as the unintegrated gluon distribution, which evolves according to the Balitsky-Fadin-Kuraev-Lipatov (BFKL) evolution equation [27][28][29]. By contrast, the TMD factorization formula emerges from (2.4) when k t ∼ Q s (x 2 ), it is formally obtained by replacing H ag * →cd (P t , 0); in that regime, (leading-twist) non-linear effects are important, and induce significant differences between the gluon TMDs.
Starting from the TMD formula, restoring the off-shellness of the small-x gluons in the hard factors and hereby obtaining the ITMD formula is equivalent to performing an all-order resummation of power corrections in |k t |/|P t |, known as kinematical-twists corrections [18]. Furthermore, the difference between the ITMD formula and the more complete CGC formulation represents corrections of the genuine-twists kind [19], that should become important when P t ∼ Q s (x 2 ). Diagrammatically, those genuine-twist corrections come from Fig. 1 (b), meaning 3-body and 4-body terms after squaring. At the cross-section level, all contributions in Fig. 1 involve 4 Wilson lines (4 fundamental ones in the case of the qq final state considered here), but the 3-(resp. 4-) body contribution involves 3 (resp. 4) derivatives and 3 (resp. 4) different transverse positions, while the ITMD cross-section is a two-body contribution which involve 2 derivatives and 2 different transverse positions, as is explicit in (2.8).

CGC framework for forward qq pair production
In this subsection, we recall the CGC formalism for qq pair production in dilute-dense collisions. In the amplitude and complex conjugate amplitude, the incoming gluon from the dilute projectile may split into the qq pair before or after the interaction with the dense target, as pictured in Fig. 2. Fundamental Wilson lines describe the interaction for quarks, and adjoint Wilson lines for gluons. As a result, the cross-section involves four contributions: a correlator of four fundamental Wilson lines, S (4) , corresponding to interactions happening after the gluon splitting into the qq pair, both in the amplitude and the complex conjugate amplitude; a correlator of two adjoint Wilson lines, S (2) , corresponding to interactions taking place before the gluon splitting, both in the amplitude and the complex conjugate amplitude; two correlators of three Wilson lines, S (3) , for the interference terms. Denoting p the momentum of the incoming gluon, the cross-section reads [24]: denote the transverse positions of the final-state quark in the amplitude and the conjugate amplitude, respectively, and denote the transverse positions of the final-state antiquark in the amplitude and the conjugate amplitude, respectively. u −u is conjugate to the hard momentum P t = (1−z)p 1t −zp 2t , and v − v is conjugate to the total transverse momentum of the produced particles k t = p 1t + p 2t . The S (i) Wilson line correlators are given by (2.16) The functions ϕ λ αβ denote the g → qq splitting wave functions. In the limit of massless quarks, the wave function overlap is simply given by The three scales Q s , |k t |, and |P t | are characterizing the kinematics for the dijet production. It is instructive to consider the two limits Q s |k t |, |P t | and Q s , |k t | |P t | in the CGC framework.
It was shown in Ref. [16] that in the Q s k t ∼ P t limit, the formula (2.11) reduces to It corresponds to the BFKL limit of the CGC and is referred to as the HEF formula. It has been extensively studied in the literature [30][31][32][33][34] (where the gluon TMD is denoted by F g/A = πΦ g/A due to a different normalization convention) 2 . Its domain of validity corresponds to jets produced away from the back-to-back region, where the small-x 2 gluon is hard, and saturation effects are negligible. However, provided that we are dealing with forward jets, linear small-x effects are still relevant [34].
In the meantime, it was shown in Ref. [24,25] that in the Q s ∼ k t P t limit, the formula (2.11) becomes This is a TMD factorization, obtained from the CGC by extracting the leading 1/P t power. Its domain of validity corresponds to nearly back-to-back jets production. In our small-x context (forward jets), saturation effects in Q s /k t must be accounted for here, without them the TMDs all coincide. We note that the TMD approach has been previously extensively studied in the literature [22,23,[36][37][38][39][40][41] in a broader context than small-x physics, in which case the process dependence of the TMDs is simply non-perturbative (Qs is not large enough compared to Λ QCD ). The ITMD factorization formula (2.4) (as well as those for the other two channels) was built in order to contain both those expressions as its limiting cases, as therefore be valid regardless of the magnitude of |k t |. In the first case, it is so because in the Q s In the second case, it occurs because in the |k t | |P t | limit the coefficient in front of F adj becomes −2z(1 − z). We note that, any systematic improvements of the HEF or TMD factorization frameworks in perturbation theory, which may be obtained in the future, could be implemented in the ITMD factorization formula as well.
Finally, the difference between the CGC formula (2.11) and the ITMD formula (2.6) was clarified recently [18]. The CGC amplitude pictured in Fig. 2, whose square leads to (2.11), can be rewritten in an alternative way using an expansion in the dipole sizes conjugate to P t , which corresponds to a twist expansion. The leading contribution represented in Fig. 1 (a) (whose square leads to the ITMD formula (2.6)), is made of the leading 1/P t term (the TMD term extracted in [24,25]) and an all-order resummation of a sub-set of higher-order terms, the so-called kinematical twists of order O(k t /P t ) (whose effect at the cross-section level is to restore the off-shellness of the gluon in the hard factor while leaving the leading-twist TMD structure unchanged). The remaining higher-order contributions represented in Fig. 1 (b), of order O(Q s /P t ), represent the difference between the CGC and the ITMD formula, they are known as genuine twists terms. Our goal now is to estimate the magnitude of that difference. To do that, we shall consider the large-N c limit, in which case the CGC framework becomes tractable, especially with the gg → qq channel.

ITMD/CGC comparison in the large-N c limit
To enable an ITMD/CGC comparison easily, let us simplify Eqs. (2.6) and (2.11). As for the multi-point correlators in (2.11), in terms of the fundamental Wilson lines, we can write down those as Indeed, Eqs. (2.21), (2.22), and (2.23) are yet too complicated to handle analytically and numerically. Therefore, we shall utilize the so-called Gaussian approximation of the CGC [42][43][44][45][46][47][48]. The essential point is to assume that all the color charge correlations in the target stay Gaussian throughout the evolution. On top of that, for simplicity, we shall work in the large-N c limit. In addition to dropping the explicitly large-N c suppressed terms, this allows to write a correlator of a product of traces as the product of single trace correlators.
Thus, the combination inside the brackets · in Eq. (2.11) can be cast into 25) in terms of only the two-point function (dipole amplitude) S qq (x, y; Then, introducing the dipole amplitude in the momentum space, and neglecting the b dependence of F for simplicity, the second and third lines of Eq. (2.11) simplify into where S ⊥ represents the transverse area of the target. Finally, with these approximations the cross section for producing a pair of q at y 1 with p 1t andq at y 2 with p 2t in the forward rapidity region is given by In the massless quarks limit, this is simply given by where we have used the identity: (2.30) Therefore, provided the large-N c limit, the CGC formula for forward dijet production reads We have performed the change of variable q t → k t − q t and then wrote (1 − z)k t − q t = (1 − z)(k t − q t ) − zq t before squaring (this choice for writing (2.31) will make for easier comparisons with the ITMD formula). Also, we have put N c /(2C F ) → 1 in the overall prefactor.
Next, let us examine the ITMD framework by using the same approximations as illustrated above. With our approximations, the ITMD framework for forward dijet production now involves only two gluon TMDs, and from (2.8), they can be written as: The forward dijet cross section in the ITMD framework is then given by (2.34) To reach the second line of Eq. (2.34), we have used the change of variable q t → k t − q t . This can now be compared with Eq. (2.31).
Let us emphasize the purpose of this paper again. In this section, we have highlighted the difference between the ITMD and CGC frameworks analytically, using the Gaussian truncation and the large-N c limit: the hard scattering part in Eq. (2.34) differs from the one in Eq. (2.31). Our interest now is to estimate the genuine twist corrections absent in the former but present in the latter. In the following section, we shall further examine that, numerically.
Before, it is worthwhile to give the HEF and TMD limits using the same simplifications, as we shall numerically evaluate them later as well. In the HEF limit where |k t | Q s ,  Meanwhile, in the TMD limit, from (2.20) one obtains dσ(pA → qqX) dy 1 dy 2 d 2 p 1t d 2 p 2t TMD = α s S ⊥ 2π 2 z(1 − z)P qg (z) (2.36)

Numerical setup and results
In order to exemplify the accuracy of the ITMD framework, we evaluate the azimuthal angle correlation in forward qq dijet production with the ITMD formula (2.34) and with the CGC formula (2.31) in p + p and p + A collisions and compare these results.

Setup
Let us first specify the setup for our numerical calculations. We assume that the prefactors α s S ⊥ in the formulas (2.34) and (2.31) are common constants and cancel out when we take a ratio of these dijet cross-sections. For the collinear gluon distribution f g/p on the projectile side, we use the parametrization CTEQ6M [49] with the factorization scale set to µ = (|p 1t | + |p 2t |)/2.
For the small-x gluons F (x 2 , k t ) on the dense target side, we include the x-evolution effects by adopting a numerical solution to the BK equation [10,11]: where Y ≡ ln(1/x 2 ) is the evolution rapidity, r ⊥ = r 1⊥ + r 2⊥ the size of a parent dipole.
In the Gaussian truncation, S qq = (S BK ) 1−1/N 2 c [43], therefore in the large-N c we simply use S BK to obtain the Fourier transform F using (2.26). The possible impact parameter dependence of the dipole amplitude is neglected here. We employ the kernel K with running coupling corrections, which was derived in Ref. [50], and we adopt the one-loop running coupling constant in coordinate space α s (r 2 ⊥ ) = 9 4π ln 4C 2 The parameter a is a smooth cutoff to make the coupling finite in the large-dipole limit: α s (r ⊥ → ∞) = 0.5 [51,52]. Our result on dijet production here is insensitive to this particular choice. For our purpose of ITMD/CGC comparison, we take as the initial condition of the BK equation the McLerran-Venugopalan (MV) type model [53,54] of the form: where x = x 0 denotes the start of the small-x evolution, which we take equal to 0.01. Other parameters are set as Q 2 0,p = 0.2 GeV 2 and Λ = 0.241 GeV for the proton target, as indicated by global fitting analysis of deep inelastic scattering small-x data with the BK equation [55][56][57]. Figure 3 displays the gluon TMDs, F gg (x 2 , k t ) and F adj (x 2 , k t ), obtained by solving the BK equation with the MV initial condition. The BK evolution contains two competitive effects, the gluon branching and merging, which results in the increase (decrease) of the gluon distributions in high (low) k t region with decreasing x 2 . In the following, we define the saturation scale Q s (x 2 ) by the peak position of the gluon TMD, F adj (x 2 , k t ), as a function of k t for fixed x 2 , which is indicated with vertical dash-dotted arrows in Fig. 3 (hence Q s (x 0 ) 0.5 GeV is slightly different from Q 0 ). The Q s (x 2 ) value increases as x 2 decreases from x 2 = 10 −2 , 10 −4 to 10 −6 . Those results are consistent with previous studies [17,58].
For a heavy nuclear target, we replace the initial Q 0 value at x = x 0 by where we have introduced a parameter c [33]. In Ref. [59] it is shown that c ≈ 0.25 − 0.5 yields a reasonable fit to the nuclear structure function F 2,A (x, Q 2 ) at x = 0.0125 measured by New Muon Collaboration. Indeed, the CGC model calculation with a smaller value ofĉ ∼ 3 (c ∼ 0.5) has resulted in more reasonable description of forward heavy-flavor production as well as quarkonium production in p + A collisions [52,[60][61][62][63] compared to the early predictions withĉ = 4-6 [51,64]. In this paper, we chooseĉ in the range of 2 ≤ĉ ≤ 3 for the initial saturation scale in heavy nuclei, Pb (A = 208) and Au (A = 197).

Kinematics
The total and relative momenta squared, (2.3) and (2.5), of the quark at y 1 with p 1t and the antiquark at y 2 with p 2t read, respectively,

4)
|P t | 2 = |p 2t | 2 |p 1t | 2 (|p 1t |e y 1 + |p 2t |e y 2 ) 2 e 2y 1 + e 2y 2 − 2e y 1 +y 2 cos φ , (3.5) where φ is the azimuthal angle between p 1t and p 2t . For definiteness, we set |p 1t | = |p 2t | = |p t | and then By changing the azimuthal angle φ, we can scan k t values from Q s ∼ |k t | |p t | to Q s |k t | ∼ |p t |. The relative momentum P t depends on the rapidities, y 1,2 . For the equal rapidity, y 1 = y 2 = y > 0, it simplifies to and for the rapidities with a large gap, y 2 y 1 > 0, it is approximated by which is almost independent of φ as e −(y 2 −y 1 ) 1. Figure 4 shows |k t | (dotted) and |P t | (dashed) as a function of the azimuthal angle φ with fixed (a) |p t | = 10 GeV and (b) 40 GeV for the pair at a common rapidity, y = 3, at √ s = 7 TeV. In the lower panels of Fig. 4 shown are the same plots but with rapidity difference, (c) y 1 = 1 and y 2 = 3, and (d) y 1 = 1 and y 2 = 5, with fixed |p t | = 10 GeV. The x 2 value is fixed by Eq. (2.2), and the corresponding saturation scale Q s (x 2 ) is indicated with dash-dotted line in each plot. The HEF formula is justified for Q s |k t | ∼ |P t |, i.e., away from the correlation limit. On the other hand, the TMD factorization formula applies to the kinematical region, |k t | ∼ Q s |P t |. We remark here that the ITMD formula will be less accurate due to genuine highertwist corrections in powers of Q s /|p t | when the separation of the scales |p t | and Q s becomes marginal by lowering |p t |, while the CGC formula is valid as long as Q s (x 2 ) Λ QCD . Before closing this subsection, we comment on the singularity appearing in the integrand at q t − p 2t = 0 in the CGC formula (2.31). It is not present in the ITMD formula and entirely pertains to the genuine-twists terms of Fig. 1 (b). It corresponds to the initial collinear gluon splitting into collinear quark/antiquark, which then independently pick up their transverse momentum from the two gluons (one with p 1 and the other with p 2 ), which indeed requires a two-body contribution at the amplitude level. In principle, this logarithmic divergence should be absorbed into a double-parton-distribution contribution not considered here [65]. For simplicity however, in this work we regularize it by adding a small mass term in the numerator as 1/((q t − p 2t ) 2 + m 2 ) in Eq.  1/p 2 2t with 1/(p 2 1t + m 2 ) and 1/(p 2 2t + m 2 ) in Eq. (2.34) for consistency. We examined the m-dependence of our numerical results by comparing results of m = 1 MeV Λ QCD and 100 MeV ∼ Λ QCD . We found no significant change in the ratios of the dijet cross-sections of the ITMD and CGC formulas at the LHC energy when p ⊥ ∼ 30 or 40 GeV. However, the change becomes noticeable around φ ∼ 0 at lower p ⊥ in both RHIC and LHC energies. In the following calculations, we will take m = 100 MeV and study the region φ > π/2.

ITMD/CGC ratio in p + p
Our focus is on azimuthal angle correlation in forward quark dijet production. We will compute the dijet yield dN (pp/pA → qqX) dy 1 dy 2 dp 1t dp 2t dφ ≡ 2πp 1t p 2t S ⊥ dσ(pp/pA → qqX) dy 1 dy 2 d 2 p 1t d 2 p 2t . The cross-section depends on the relative angle φ, not on individual angles of p 1t,2t due to the rotational symmetry of the dijet production. Note that the qq dijet yield is given by N = σ qq /S ⊥ with the assumption S ⊥ ≈ σ inel between the effective transverse area and the inelastic cross section. Both the CGC and ITMD formulae contain the TMD and HEF limits within the appropriate kinematics, the difference between them represents genuine higher-twist contributions, present in the CGC results but absent in the ITMD case, where only kinematical-twist contributions in |k t |/|P t | are resumed. In terms of the CGC formula, that difference comes from power corrections in the dipole size expansion 3 .
In order to quantify the genuine higher-twist effects, we compare results of the ITMD formula and of the CGC by taking the ratio of the former to the latter. In Fig. 5, we show the ratio R as a function of φ for the pair of the common rapidity y 1 = y 2 = y at (a) |p t | = 40 GeV and (b) 10 GeV. The ITMD/CGC ratio R for y = 1 (thick black line) in Fig. 5 (a) is consistent with unity over the whole range of φ studied here, and for y = 4 it lies barely below unity as Q s (x 2 ) becomes larger. Other ratios of TMD/CGC (blue dashed) and HEF/CGC (red dashed) deviate from unity outside of their respective domain of applicability, as is expected. Indeed, the dijet production cross-section in the HEF formula unphysically vanishes dσ HEF (k t → 0) → 0 in the back-to-back limit (see Eq. (2.35)). On the other hand, the TMD formula, which ignores the off-shellness of the partons in the hard matrix factors, underestimates the cross-section for φ away from the back-to-back region.
At the lower |p t | = 10 GeV (Fig. 5 (b)), the ITMD/CGC ratio R shifts below unity by 5 to 15% around φ ∼ π/2 as the rapidity y increases from y = 1 to 4. This deviation can be understood as a result of the increase of the power corrections mentioned earlier. We stress here that the ITMD formula approximates the CGC result uniformly over the region of φ with 5-15% accuracy. The genuine higher-twist corrections become more important outside the back-to-back region, while it is negligible around the back-to-back limit φ = π. We will investigate the power corrections further in the last subsection below.
In Fig. 6, we study the energy dependence of the ITMD/CGC ratio R by showing the cases of √ s = 2.76 TeV (thick) and 13 TeV (thin) for (a) |p t | = 40 GeV and (b) 10 GeV. For the larger p t (a), we find that the ITMD approximation is very accurate there, almost the same result as in Fig. 5, which indicates that the corrections in Q s /p t are well suppressed there. For the lower p t (b), the ITMD/CGC ratio deviates from unity and the depletion becomes more significant with increasing collision energy √ s = 2.76 to 13 TeV (i.e., increasing We also examine dijet production with a rapidity separation in the cases, (y 1 , y 2 ) = (1, 3) and (1,5), as shown in thick and thin lines, respectively, in Fig. 7 for (a) |p t | = 40 GeV and (b) 10 GeV. From Fig. 7 (a), we see that the ITMD formula well approximates the CGC result, and interpolates the results of TMD and HEF formulas uniformly over the range of π/2 φ π. The TMD estimate becomes accurate near the back-to-back region φ ∼ π, while it is less accurate for φ away from it. For larger y 2 , this approximation becomes better. On the other hand, unsurprisingly, the HEF formula fails to reproduce the CGC result in a wider region of φ around ∼ π irrespective of the y 2 value.
At the lower |p t | = 10 GeV (Fig. 7 (b)), the ITMD formula interpolates still smoothly from the TMD to the HEF results with decreasing φ from the back-to-back limit φ ∼ π. However, the value of the ITMD/CGC ratio lies significantly below unity in the non-backto-back region, which reflects the size of the genuine twist corrections.

ITMD/CGC ratio in p + A and nuclear modification factor
The saturation scale Q 2 sA in a heavy nucleus will be enhanced by a factor ofĉ ∝ A 1/3 compared to Q 2 sp , as discussed in Sec. 3.1, and therefore it is valuable to analyze the nuclear dependence of the ITMD/CGC ratio in forward dijet production in p + A collisions. We plot in Fig. 8 the ratios R in p + p (ĉ = 1) and p + A (ĉ = 2.5) collisions at √ s = 5.02 TeV for |p t | = 40 and 20 GeV ((a) and (b)), and for |p t | = 10 and 5 GeV ((c) and (d)). From the comparison of the R ratios in p + p and p + A collisons at |p t | = 40 and 20 GeV in Fig. 8 (a) (b), we find that the deviation of the ratio R from unity becomes more noticeable in p + A collisions and for the lower |p t | = 20 GeV, which indicates the enhanced power corrections of Q s /|p t | at lower p t . At yet lower values |p t | = 10, 5 GeV, the deviation becomes significant even in p + p case (Fig. 8 (c)), and is more profound in p + A case (Fig. 8  (d)). In these cases, the ITMD is no longer a good approximation to the CGC. Figure 9 shows the results at the RHIC energy, √ s = 200 GeV. Since the dijet production formulas, Eqs. (2.34) and (2.31) premise that x 2 is small, x 2 ≤ x 0 = 0.01, the jet momentum |p t | is accordingly limited to the lower values, and here we take y = 2 and |p t | = 5 and 3 GeV. Although Q s becomes smaller at RHIC, the ratios R for |p t | = 5 and 3 GeV at √ s = 200 GeV in Fig. 9 deviate from unity similarly to those for |p t | = 10 and 5 GeV at √ s = 5.02 TeV in Figs. 8 (c) (d). This result shows the dijet production at the RHIC energy is sensitive to genuine higher-twist corrections in Q s /|p t |, which are not included in the ITMD formula. Next let us discuss the so-called nuclear modification factor R pA , for forward dijet production (dP.S. = dy 1 dy 2 dp 1t dp 2t dφ):  where S p,A ⊥ denote the effective transverse areas of the proton and nucleus targets, respectively. A p + A collision should be presumably regarded as a superposition of p + p collisions for the high momentum limit |p t | → ∞ (at φ = π), and then R pA → 1 is expected. To assure this constraint, we normalize the effective transverse area in our model calculations as A modification of R pA from unity signals the presence of nuclear effects. Figure 10 demonstrates R pA of the quark dijet production at y = 3 at √ s = 5.02 TeV for jet momentum |p t | = 40, 20, 10, and 5 GeV. Colored bands depict the uncertainty estimated by the change of the results when the initial saturation scale of the nucleus Q 2 s0,A =ĉ Q 2 s0,p is varied in the range ofĉ = 2-3. At |p t | = 40 GeV ( Fig. 10 (a)), the CGC (blue solid), and ITMD (red dashed) formulas give the same prediction for R pA . The prediction is consistent with unity over a wide range of φ and is suppressed only in the TMD regime in the vicinity of φ = π, where the total transverse momentum of the dijet becomes small and comparable to the saturation scale: |k t | Q s (x 2 ). The intrinsic transverse momentum of the gluons, which is of the order of Q s , is larger in the heavy nucleus than in the proton and smears the azimuthal angle correlation. The region of the suppression should be characterized by δφ = |π − φ| Q sA /|p t | which we indicate with a vertical dash-dotted arrow in Fig. 10 (a). Decreasing the jet momentum |p t | from (a) 40 GeV down to (d) 5 GeV in Fig. 10, we find that the suppression of R pA appears in a wider range of φ on the away side. This is because, for lower |p t | at fixed y, the relevant x 2 is smaller and Q s (x 2 ) is larger accordingly, and therefore the region of δφ Q s /|p t | gets wider.
We also notice that the difference between the ITMD and CGC results increases as |p t | decreases. In the suppression regime, near φ = π, the differences stay rather small, but at moderate φ away from φ ∼ π, those differences can get large and in fact, in that regime the CGC formula exceeds unity while the ITMD one stays suppressed. This qualitative change is caused by the genuine higher-twist corrections ( Fig. 1 (b)) present in the CGC formula. They contribute significantly to quark dijet production at moderate values of φ (where the ITMD cross-section is not very large), and contribute even more so in p + A collisions compared to p + p collisions, due to the bigger saturation scale in the former case. That creates an enhancement of R pA .
At the RHIC energy √ s = 200 GeV, a similar suppression of R pA is seen in the back-toback region around |δφ| < Q sA /|p t | in Fig. 11, reflecting the larger intrinsic k t in the nuclear target. In contrast, at φ away from π, both the CGC and ITMD results show enhancements of R pA . The larger discrepancy between the ITMD and CGC results in Fig. 11 (b) than in (a) is again a manifestation of the larger genuine twist corrections to dijet production at |p t | = 3 GeV compared to |p t | = 5 GeV. The enhancement of the ITMD cross section at φ away from π is not surprising actually; at RHIC energies we are simply sensitive to our initial conditions: if we plot the ratio of the gluon TMD F gg (x 2 , k t ) of the heavy nucleus to that of the proton at x 2 ∼ x 0 , it does show a Cronin-like peak structure as a function of k t . See Fig. 12. Indeed, Fig. 12 compares R pA of the quark dijet production cross-section obtained with the CGC and ITMD formulas, together with R pA of the gluon TMDs F gg and F adj (in order to highlight the higher-twist effects, we choose the lower values of the jet momentum |p t | = 10 and 5 GeV, resp. at the LHC (a) and RHIC (b) energies). We see that the cross-section ratio obtained with the ITMD formula is roughly proportional to the ratio of F gg (and, away from φ = π, to F adj also, when that TMD is no more proportional to k 2 t , see Fig. 3). The contribution of multi-body scattering diagrams in the CGC to higher-twist corrections was also addressed in Refs. [42,66], as k t -factorization breaking effect for quark- antiquark pair production. In that analysis, what was studied was the difference between the CGC and HEF formulae, which contains two types of HEF (or k t ) factorization breaking contributions: leading-twist saturation corrections in Q s /|k t | and genuine-twist saturation corrections in Q s /|p t |. In the present work, by employing the ITMD framework, we are now able to include the former in the baseline, and isolate the latter as the difference between the CGC and ITMD formulae. To illustrate our findings, Fig. 13 displays the dijet production yield at the LHC and RHIC as a function of |k t |/|P t |. We find that at small values of |k t | (around |k t |/|P t | = O(0.1) or smaller), the leading-twist saturation corrections are responsible for the (rather large) difference between the CGC and HEF cruves, as the genuine-twist corrections are negligible (since the ITMD and CGC curves coincide). By contrast, the genuine-twist saturation corrections become visible when |k t |/|P t | 1, where the HEF and ITMD cross-sections are equal (implying negligible leading-twist saturation corrections), but both different from the CGC one. The figure shows the maximal size of the genuine higher-twist effects, which become more visible with heavy nuclear target (largê c).

Summary
We have compared quantitatively in detail the result of the ITMD formula to that of the CGC formula for forward qq dijet production in p + p and p + A collisions. We assumed that the typical transverse momentum of a hard jet P t is much bigger than the saturation scale of the target Q s , but considered arbitrary values of k t , the transverse momentum imbalance of the quark-antiquark pair. First, Sec. 2 has recaptured the differences and similarities between the two frameworks in describing the forward dijet production cross-section. The ITMD formula (2.6) contains three kinds of leading-twist small-x gluon TMDs, but two of them, F gg and F adj , are relevant in the large-N c approximation. At small k t (the TMD regime), the differences between those distributions, see Fig. 3, is the result of an all-order resummation of saturation corrections in Q s /|k t |, while the hard factors incorporate an all-order resummation of kinematical twists in |k t |/|P t |, resulting in a proper matching to BFKL at large k t (the HEF or k t -factorization regime). The CGC formula (2.11) involves 2-, 3-and 4-point correlators of Wilson lines; it contains the full ITMD formula and on top resums the genuine higher-twist contributions in Q s /|P t |.
Using the Gaussian truncation, however, one can make the two formulae look rather similar: both involve convolutions of the qq dipole amplitude in momentum space (Eq. (2.26)) with itself and with hard parts. In the ITMD case (2.34), those convolutions are simply the gluon TMDs (2.33). In the CGC case they are more involved (2.31) as they include the genuine higher-twist contributions. Those come from multi-body correlators (e.g. (2.21)), but our approximations have allowed us to write them in terms of the function F . The genuine high-twists are suppressed in high-p t dijet production, i.e., |P t | Q s , in which case the ITMD formula represents a good approximation to the CGC framework.
In Sec. 3, we have demonstrated the quantitative difference between the two formulas for forward quark dijet production by evaluating the azimuthal dijet correlation in p + p and p + A collisions at collider energies. We have confirmed that the ITMD formula, which interpolates between the TMD and HEF formula, gives the same prediction as the CGC for the dijets with p t ∼ 40 GeV at the LHC energy, where the higher-twist genuine corrections are suppressed. As p t is decreased, some difference is seen and amount to around 5-15% for p t ∼ 10 GeV at moderate φ away from the back-to-back limit. We can regard that amount as the highest estimation of the genuine twist effect for the qq dijet correlation, as well as for the other dijet channels for which those estimations would be more involved.
The nuclear modification factor R pA in p + A collisions shows a dip structure around the back-to-back region of φ in both the frameworks, resulting from leading-twist saturation effects in nuclear versus proton targets, and reflecting the intrinsic k t of the gluons, which is of the order of Q s . For p t 10 GeV at moderate φ away from the back-to-back limit at the LHC, the ITMD gives a suppression while the CGC formula yields an enhancement. We attribute this difference to a nuclear enhancement of the genuine-twist contributions, i.e., the higher-body multiple scattering effects included in the CGC formula. The effects are more substantial at lower p t and with the denser nuclear target than higher p t with the dilute one.
When the ITMD formula is used to evaluate the forward dijet production cross-section at moderate p t for the study of gluon saturation, one should be aware of the fact that this framework lacks those genuine twist effects. We note that the studies which are restricted to the TMD regime near φ = π, e.g. to isolate the contribution of polarized gluons (relevant when massive quarks are considered [26,67], for dijets in deep inelastic scattering [68][69][70][71][72] or for three-particle production [73][74][75]) or to implement a Sudakov resummation [76,77], are rather safe provided the p t 's are not too low.
It would also be interesting to examine whether those effects are experimentally measurable, provided that the gluon TMDs could be determined with good accuracy in other processes. For this purpose, we need to take account of the effects of jet fragmentation and also other effects in jet identification algorithm and efficiency cuts, and so forth. We leave those as future work.