Forward di-jet production in p+Pb collisions in the small-x improved TMD factorization framework

We study the production of forward di-jets in proton-lead and proton-proton collisions at the Large Hadron Collider. Such configurations, with both jets produced in the forward direction, impose a dilute-dense asymmetry which allows to probe the gluon density of the lead or proton target at small longitudinal momentum fractions. Even though the jet momenta are always much bigger than the saturation scale of the target, $Q_s$, the transverse momentum imbalance of the di-jet system may be either also much larger than $Q_s$, or of the order $Q_s$, implying that the small-$x$ QCD dynamics involved is either linear or non-linear, respectively. The small-$x$ improved TMD factorization framework deals with both situation in the same formalism. In the latter case, which corresponds to nearly back-to-back jets, we find that saturation effects induce a significant suppression of the forward di-jet azimuthal correlations in proton-lead versus proton-proton collisions.


Introduction
Measurements of forward particle production in high-energy hadronic collisions provide unique opportunities to study the QCD dynamics of the non-linear parton saturation regime [1]. Such processes, in which, for kinematical reasons, high-momentum partons from one of the colliding hadrons mainly scatter with small-momentum partons from the other, are called dilute-dense collisions. Indeed, the density of the large-x partons in the projectile hadron is small, while the density of the small-x gluons in the target hadron is large, and the former, well understood in perturbative QCD, can be used to probe the dynamics of the latter. This is true already in proton-proton collisions, although using a target nucleus does enhance the dilute-dense asymmetry of such collisions.
RHIC measurements have provided some evidence for the presence of saturation effects in the data, the most compelling of which is the successful description of forward di-hadron production [2][3][4], using the most up-to-date theoretical tools available at the time in the Color Glass Condensate (CGC) framework [5,6]. In particular, this approach predicted the suppression of azimuthal correlations in d+Au collisions compared to p+p collisions [7], which was observed later experimentally [8,9].
The CGC effective theory provides a tool to compute observables when non-linear QCD dynamics must be taken into account. It effectively describes, in terms of strong classical fields, the dense parton content of a hadronic/nuclear wave function, at small longitudinal momentum fraction x. The separation between the linear and non-linear regimes is characterized by a momentum scale Q s (x), called the saturation scale, which increases as x decreases, and roughly scales as A 1/6 . The CGC description of dilute-dense collisions from first principles is valid provided Q s Λ QCD , therefore it should work better with higher energies, as they open up the phase space towards lower values of x. In order to verify that this is the case, the CGC predictions must be extended from RHIC kinematics to the Large Hadron Collider (LHC), where the relevant observables involve high-p t jets, as opposed to individual hadrons with p t of the order of a few GeV at RHIC.
In this context, we shall consider forward di-jet production in proton-lead versus protonproton collisions. In that case, it was shown in [10] that the full complexity of the CGC machinery is not needed. Indeed, for the di-hadron process at RHIC energies, no particular ordering of the momentum scales involved is assumed in CGC calculations, while at the LHC one can take advantage of the presence of final-state partons with transverse momenta much larger than the saturation scale to obtain simplifications. On the flip side, different complications -left for future studies -are expected to arise due to QCD dynamics relevant at large transverse momenta and not part of the CGC framework, such as Sudakov logarithms [11][12][13][14] or coherence in the QCD evolution of the gluon density [15][16][17].
There are three distinct momentum scales in the forward di-jet process. The typical jet transverse momentum P t is always one of the hardest scales, and it is much bigger than the saturation scale Q s , which is always one of the softest scales. The third momentum scale is the total transverse momentum of jet pair k t , which also corresponds to the transverse momentum of the small-x gluons involved in the hard scattering. Depending on where k t sits with respect to P t and Q s , the full CGC formulation simplifies either to the high energy factorization (HEF) framework [18,19] or to the (small-x limit of the) transverse momentum dependent (TMD) factorization framework [20].
The HEF framework is recovered from the CGC when Q s k t ∼ P t [10]. In that case, non-linear effects are absent, and the description of forward di-jets involves off-shell hard matrix elements, along with a single TMD gluon distribution for the small-x target (also called unintegrated gluon distribution in the literature). The TMD framework is recovered from the CGC when k t ∼ Q s P t [21][22][23]. In that case, some non-linear effects do survive, and the description of forward di-jets involves several TMD gluon distributions, each associated to a sub-set of the hard matrix elements, but those are on-shell. In Ref. [10], we proposed an interpolating formula between those two limits, applicable for P t Q s regardless of the magnitude of k t , which is more amenable to phenomenological implementations than the CGC expression which, as we pointed out, also contains both the HEF and TMD limits. Not unexpectedly, it involves several unintegrated gluons distributions each associated to a sub-set of off-shell matrix elements. The goal of this paper is to provide a numerical implementation of that new formulation, dubbed improved TMD (ITMD) factorization.
The off-shell matrix elements needed to compute the forward di-jet process have all been calculated in [10], but evaluating all the necessary gluon TMDs is not straightforward. Very recently, they have been obtained from a numerical simulation of the non-linear QCD evolution in the leading ln(1/x) approximation [23], that is from the Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) [24][25][26][27][28] equation. However, further work is required before those TMDs can be incorporated into a cross section calculation. Therefore, in the present work, we shall stick to a mean-field type approach in which all the gluon distributions needed can be related to each other, and obtained from the simpler Balitsky-Kovchegov (BK) equation [29,30]. A detailed comparative study using solutions of the different extensions of the original BK equation is left for future work. The version that we shall use in this work is known as the KS gluon distribution [31]. It incorporates the running of the QCD coupling, non-singular pieces (at low x) of the DGLAP splitting function, a sea-quark contribution, and resums dominant corrections from higher orders via a kinematic constraint [32,33].
By comparing the forward di-jet production cross sections in proton-lead and protonproton collisions, we can clearly see the onset of parton saturation effects, as we go from a kinematical regime in which k t ∼ P t towards one where k t ∼ Q s , and we obtain a good estimation of the size of those effects where they are the biggest, which is for nearly back-toback jets. We note that probing non-linear effects of similar strength with single-inclusive observables requires to make the only transverse momentum involved in those processes of the order of the saturation scale, which may not be easy experimentally. With di-jets, assuming P t ∼ 20 GeV and k t ∼ Q s ∼ 2 GeV, we can reach R pP b ∼ 0.5.
The paper is organized as follows. In section 2, we recall the essence and the ingredients of the ITMD factorization formula for forward di-jets in dilute-dense collisions. In section 3, we introduce the mean-field approximation that allows us to express the various gluon TMDs in terms of the solution of the BK equation. In section IV, we present numerical results for the proton and lead gluon TMDs obtained with the KS gluons, and compared them with analytical expressions obtained in the GBW model. In section V, we present our results for forward di-jet production in p+p and p+Pb collisions at the LHC, as well as nuclear modification factors R pPb . Finally, section VI is devoted to conclusions and outlook.

The ITMD factorization formula for forward di-jets in dilute-dense collisions
We consider the process of inclusive forward di-jet production in hadronic collisions where the four-momenta of the projectile and the target are massless and purely longitudinal. The longitudinal momentum fractions of the incoming parton from the projectile, x 1 , and the gluon from the target, x 2 , can be expressed in terms of the rapidities (y 1 , y 2 ) and transverse momenta (p t1 , p t2 ) of the produced jets as (2.2) 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 qg → qg , gg → qq , gg → gg .   ag→cd in the large-N c limit. The finite N c expressions can be found in [10]. the dense target are described by TMD distributions Φ g/A (x 2 , k t ) . Indeed, the momentum of the incoming gluon from the target is not only longitudinal but also has a non-zero transverse component of magnitude which leads to imbalance of transverse momentum of the produced jets: where P t is the hard scale of the process, related to the individual jet momenta P t ∼ |p 1t |, |p 2t |. By contrast, the value of k t can be arbitrary. The ITMD factorization formula reads [10] dσ pA→dijets+X It involves several gluon TMDs Φ (i) ag→cd (2 per channel), with different operator definitions, that are accompanied by different hard factors K (i) ag * →cd . Those where computed in [10] using either Feynman diagram techniques, or color-ordered amplitude methods, and they are given in Table 1 in terms of the Mandelstam variables of the 2 → 2 parton level process. They encompass the improvement over the TMD factorization formula derived in Ref. [22] where the matrix elements were on-shell and a function of P t only. The gluon TMDs are normalized such that and their precise operator definitions can be found in [10]. As emphasized in the introduction, formula 2.6 coincides with CGC expressions in two important limits. They both reduce to the TMD factorization formula when Q s ∼ k t P t and to the HEF formula when Q s k t ∼ P t : • The TMD factorization formula with k t dependent gluon distributions and on-shell matrix elements is simply obtained form 2.6 after simplifying K The derivation of this expression from the CGC framework was done in [22] the in large-N c limit, and in [23] for the finite N c case. However the TMD approach had been previously extensively studied in the literature [20,[34][35][36][37][38][39][40], and in a broader context than small-x physics.
• Obtaining the HEF formula with a single gluon TMD and off-shell matrix elements from Eq. 2.6 relies on the fact that up to power corrections, all the gluon TMDs coincide in the large k t limit: Then, denoting the HEF formula is This expression is also obtained from the CGC framework in the dilute target limit [10], and has also been extensively studied in the literature [13,19,31,41,42] (where the gluon TMD is denoted F g/A = πΦ g/A due to a different normalization convention).
We would like to point out that the ITMD factorization formula 2.6 was build in order to contain both the HEF and the TMD expressions as its limiting cases, and as such should be considered no more than an interpolating formula. We note however, that if one would be able to directly derive a factorization formula valid for Q s P t regardless of the value of k t , any additional term compared to 2.6 should vanish in both limits Q s ∼ k t P t and Q s k t ∼ P t .

The gluon TMDs in the Gaussian approximation
The goal of this paper is to provide a numerical implementation of the ITMD factorization formula, which first requires to evaluate all the gluon TMDs that enter Eq. 2.6. Let us start with the simplest of them, Φ qg→qg , also called the dipole gluon distribution and often denoted x 2 G (2) . In the small-x 2 limit, it can be related to the Fourier transform of the fundamental dipole amplitude N F (x 2 , r) where r denote the transverse size of the dipole [22,23]: The amplitude N F is defined through the CGC expectation value of the S-matrix, S F , of a quark-antiquark dipole scattering off the dense target: The dipole gluon distribution can then be written in a compact form as: and with S ⊥ denoting the transverse area of the target. In full generality, none of the other gluon TMDs can be obtained in such a straightforward manner. For instance, the Weizsäcker-Williams (WW) gluon distribution, denoted x 2 G (2) , should be obtained in the small-x 2 limit from the quadrupole operator Tr , and in general is not related to F (x 2 , k t ). Therefore, in order to simplify the evaluation of all the gluon TMDs which we need, we will resort to a mean-field type approximation.
We shall utilize the so-called Gaussian approximation of the CGC [7,[43][44][45][46][47][48]. The essence of this approximation is to assume that all the color charge correlations in the target stay Gaussian throughout the evolution: ρ(x)ρ(y) x ∝ µ 2 (x, x − y). In addition, for simplicity, we shall work in the large-N c limit. This Gaussian approximation allows to write, among other things, the WW gluon distribution in terms of an adjoint dipole: where now, S A (x, r) is an S-matrix for the scattering of a gluon dipole involving adjoint Wilson lines. The Gaussian approximation also allows to write S F = S and S A = S 2 BK , where C F and C A are the Casimirs of the fundamental and adjoint representations of SU (N c ), respectively, and with S BK denoting the solution of the BK equation. At large N c , 2 , and one can write: Then the Laplacian can be inverted as: In the large N c limit, the six gluon distributions Φ (i) ag→cd reduce to [10]: Therefore, we need an input of five gluon TMDs in our numerical calculations, the dipole gluon distribution, and four others: gg , and F (6) gg . The WW distribution is not directly one of them, but in the Gaussian approximation coupled to the large-N c limit, which ensures the factorization of CGC expectation values into single trace expectation values, those four gluon distributions can be expressed in terms of x 2 G (1) and x 2 G (2) [22]: 14)

Results for the gluon TMDs
Before we proceed with the computation of the gluon TMDs (3.2) and (3.12)-(3.16) from a solution of the BK equation, we would like to give a couple of useful and interesting results. First, we obtain the gluon TMDs in the Golec-Biernat-Wusthoff (GBW) model analytically; those results may be used for various purposes, such as checking the numerical procedure needed for the various convolution in (3.12)-(3.16). Second, we compute the high-k t behavior of the gluon TMDs in the McLerran-Venugopalan (MV), and show that it features the behavior (2.9) expected from the operator definitions of the TMDs; this was not obvious a priori, and in general not every model processes this characteristic. We recall that this behavior is necessary in order for the ITMD formula to reproduce the HEF limit when Q s k t ∼ P t . Finally, we present the gluons with we will use for our cross section calculations, obtained from the KS solution of the BK equation.

Analytical results in the GBW model
The GBW model [49] is a phenomenological model for the dipole scattering amplitude N F (x, r), that describes deep inelastic (proton) data at small-x and for moderate values of the photon virtuality. The scattering amplitude in this model is N F (x, r) = 1−exp −r 2 Q 2 s (x)/4 and the Fourier transform of S F reads: . (4.1) Using this result in formulae (3.2) and (3.8), we get x 2 G (1) and x 2 G (2) in the GBW model (see appendix A): and Note that the above result for the WW distribution can also be obtained directly from Eq. 3.4, with S A (x, r) = exp −r 2 Q 2 s (x)/2 . The expression for the WW distribution can be simplified by expressing the remaining integral in terms of the exponential integral special function, Ei(x) = ∞ x dt e −t /t: .

Large-k t behavior in the MV model
The GBW model is not a good parametrization for large transverse momenta, since not only the various TMDs do not converge to a single one at large k t , but also the exponential fall-off of each gluon TMD is unphysical. A model in which those deficiencies are corrected is the MV model [50,51]. This model comes about when the color field correlations in the Gaussian approximation are assumed to stay local µ 2 (x, x − y) → µ 2 δ(x − y). Then the saturation scale is related to the color charge density in the transverse plane of the nucleus µ 2 , integrated over the longitudinal direction: Q 2 s = g 4 C F /(2π) dz + µ 2 . In addition, the scattering amplitude in this model can be written where Λ is an infrared cut-off. The logarithmic behavior in (4.10) is only valid in the limit of small dipole sizes, but this is precisely what we need in order to study the high-k t behavior of F (x, r) and of the various gluon TMDs. As a matter of fact, the logarithm is the crucial difference between the GBW and the MV model, which restores the correct high-k t perturbative power-law behavior of the dipole gluon distribution, x 2 G (2) (x 2 , k t ), and of the WW gluon distribution, x 2 G (1) (x 2 , k t ), which both behave identically as ∼ Q 2 s /k 2 t (see for example [52] and [53]). In the appendix B, we derive the leading order term in Q 2 s /k 2 t for the remaining TMDs, by expanding Eq. (4.10) to first order in r 2 Q 2 s . We find that, expect for F (2) gg which goes to zero at leading order, they all scale the same as x 2 G (2) and x 2 G (1) : The sub-leading N c contribution to Φ (2) gg→qq (see (3.10)) is actually x 2 G (1) [10], therefore these results show that in the MV model, the behavior (2.9) is satisfied, and the ITMD formula will indeed reproduce the HEF limit when Q s k t ∼ P t . This is also true if the MV model is used as an initial condition to solve the BK equation, since the power-law fall-off (4.11) will acquire an anomalous dimension due to the small-x 2 evolution, but it will stay the same for all the gluon TMDs.

Gluon TMDs from the KS solution to the BK equation
The KS solution [31] is a solution to BK equation extended to take into account higherorder corrections relevant in order to provide realistic phenomenological predictions when somewhat large transverse momenta are involved as is the case with jets. Namely, these are corrections coming from including non-singular pieces of the gluon splitting function, kinematic constraint effects and contributions from sea quarks [33,54]. Let us already point out that the initial condition used in [31] is not the MV model, and therefore in the large k t limit, all the gluon TMDs will not exactly coincide. The mismatch is small however, probably due to the fact that the initial condition used does also effectively contain a logarithmic behavior.
The KS solution provides directly the dipole gluon TMD x 2 G (2) (x 2 , k t ), and the parameters of the initial condition are constrained by a fit to experimental data on deep inelastic scattering off protons. To deal with the nuclear case, the following formal substitution is made in the non-linear term of the equation where R A is the nuclear radius and A the mass number (A = 208 for Pb). c is a parameter that we shall vary between 0.5 and 0.75 in order to assess the uncertainty related to the strength of saturation effect in the lead nucleus compared to the proton. The nuclear dipole gluon TMD obtained in this way is also normalized to the number of nucleons A.
In order to calculate all the gluon TMDs (3.12)-(3.16) from x 2 G (2) (x 2 , k t ), we are facing the following issue. The KS solution provides directly an impact-parameter-integrated distribution, which in fact explains why the non-linear term depends on the target size. As a consequence, it is not straightforward to identify S ⊥ and obtain F (x 2 , k t ). Our procedure will be to first compute the dipole cross section σ dipole (x, r = |r|) = 2 d 2 b N F (x 2 , r) from x 2 G (2) (x 2 , k t ) by inverse Fourier transformation of Eq. (3.1), and then to define S ⊥ as its value at large r i.e. when it saturates (since in that limit N F → 1):   Figure 3: Left plot: differential cross section as a function of the azimuthal angle between the jets for p+p collisions, comparing the new ITMD approach with previously obtained HEF results. The ITMD/HEF difference, which as expected is the largest around ∆φ π is similar in p+Pb collisions, resulting in almost identical R pP b for both approaches: right plot.
We can now obtain F (x 2 , k t ) and calculate all the needed gluon TMDs. Their behavior as a function of k t is plotted in Fig. 2, both for the proton and the lead nucleus. The small mismatch between their high-k t behavior, expected due to the initial condition for the x 2 evolution, can be seen.

Numerical studies of the forward di-jet cross section
We move now to the numerical results for forward di-jet production in p+p and p+Pb collisions at the LHC. We consider a center-of-mass energy of 8. 16 TeV, and generate all our predictions with the forward region defined as the rapidity range 3.5 < y < 4.5 on one side of the detector. The two hardest jets are required to lie within this region and we also impose a cut on the minimal transverse momentum of each two jets: p t0 = 20 GeV. In such a setup, the cross section still may be divergent due to collinear singularities. These are cut-off by applying a jet algorithm on the final state momenta with a delta-phi-rapidity cut R = 0.5. Finally, we require the jets to be ordered according to increasing transverse momentum, that is we have |p t1 | > |p t2 | > p t0 .
The new factorization approach summarized in Eq. 2.6 has been implemented in two independent Monte Carlo codes avhlib [55,56] and LxJet [57]. To be more precise, the computer programs do not utilize the formula 2.6 which uses amplitudes squared and summed over polarizations and colors. Instead, the more generic color-ordered off-shell helicity amplitudes are used, as derived in [10]. This approach follows the modern direction of amplitude calculation and allows for more 'exclusive' calculations in the future (for example a study of helicity dependence or interfacing with color-flow dependent parton shower generators). The calculations of matrix elements are made keeping N c finite.
For the collinear parton distributions that enter the ITMD formula, we chose the generalpurpose CT10 set [58]. For the central value of the factorization and renormalization scale, we choose the average transverse momentum of the two leading jets, µ F = µ R = 1 2 (|p t1 | + |p t2 |). We will produce error bands corresponding to the renormalization and factorization scale uncertainties by varying the central numbers from half to twice their value.   with A = 208 for Pb. In our approach, in the absence of saturation effects, or in the case in which they are equally strong in the nucleus and in the proton, this ratio is equal to unity. If, however, the non-linear evolution plays a more important role in the case of the nucleus, the R pPb ratio will be suppressed below 1. We start by investigating the azimuthal correlations, with the azimuthal angle between the jets ∆φ defined to lie within 0 < ∆φ < π. First we compare the new ITMD approach with previously obtained HEF results in Fig. 3. For the ∆Φ distribution in p+p collisions, we see that at small angles where ideally they should match, there remains a small difference between the ITMD and HEF curves. As we anticipated, this is due to the initial condition used to obtain the KS gluons. By contrast, near ∆Φ π, we observe a large difference, as expected: the ITMD result is about a factor 3 bigger than the HEF one. The ITMD/HEF ratio is very similar in the case of p+Pb collisions, resulting in almost identical R pPb for both approaches, as also shown on the figure. For that comparison, we have parametrized the strength of the non-linear term in the evolution equation for the Pb gluon distributions (see (4.13)) with c = 0.5.
Next, we compare the ∆Φ distribution in p+p and p+Pb collisions in Fig. 4. After rescaling the p+Pb cross section by the number of nucleons, we obtain identical distributions almost everywhere. It is only for nearly back-to-back jets, around ∆φ π, that saturation effects induce a difference. This difference is better appreciated on the nuclear modification factor, which goes from unity to 0.6, as ∆φ varies from ∼ 2.7 to π. Two values of the parameter c have been considered, which makes up an uncertainty band that turns out to be rather small. This means that the uncertainty related to the value of the saturation scale of the lead nucleus does not strongly influence the predicted R pPb suppression. Finally, in Fig. 5 we display the nuclear modification factors as a function of the transverse momentum of the leading and sub-leading jet. Our conclusions are similar for these observables: the new ITMD predictions are similar to the previously obtained HEF results, due to the fact that the ITMD/HEF ratio is similar in p+p and p+Pb collisions. This means that the HEF framework, which is incorrect for nearly back-to-back jets -since in this formalism all the gluon TMDs are considered equal regardless of the kinematics -can nevertheless be safely used for R pPb calculations. The same is not true for cross section calculations. Fig. 6 shows those same nuclear modification factors but comparing the predicted suppression for two different values of the parameter c. As a function of the leading jet p t , R pPb rises up from about 0.6 for p t1 = 20 GeV to unity for p t1 = 50 GeV. However, it is interesting to note that as a function of the sub-leading jet p t , this ratio rather stays flat around 0.8.

Conclusions
In this paper, we have studied forward di-jet production in proton-proton and proton-lead collisions, using the small-x improved TMD factorization framework Eq. (2.6. We have obtained the first numerical implementation of this formalism, and the first predictions for forward di-jets at the LHC, a process which is particularly interesting from small-x point of view. Our results for the nuclear modification factors in p+Pb vs p+p collisions confirm the conclusions obtained in [41] in the HEF framework, that for nearly back-to-back jets, non negligible effects of gluon saturation are to be expected as one goes from p+p to p+Pb collisions.
This is due to the fact that in such configurations, the total transverse momentum of the jet pair is of the order of the saturation scale of the target, and even though the jet transverse momenta are individually much larger than Q s , saturation effects are not irrelevant. To obtain our predictions, we used the KS gluon distributions, and it would certainly be interesting use other extensions of the BK equation, such as for instance the rcBK gluon distribution [59], in order to compare the level of saturation effects expected.
It is important to note that so far, our results have been obtained using an impactparameter averaged nuclear saturation scale. However, the outcome of high-energy proton- nucleus collisions seems to be quite sensitive to the fact that the nucleon positions in the nucleus fluctuate event by event. We have provided predictions using two different nuclear saturation strength parameter c, but a more complete study including such nucleon-level fluctuation effects would allow to better estimate the uncertainty related to the nuclear geometry. In the meantime, our results are enough to motivate experimental measurements at the LHC. Finally, one important theoretical ingredient is still missing in our formulation: the Sudakov logarithms. Their effect should be the largest also for nearly back-to-back jets, since they are logarithms of the ratio of the hard scale to the transverse momentum imbalance of the jet pair, and therefore are expected to compete with saturation effects. The Sudakov logarithms are identical in p+p and p+Pb collisions, which means that to some extent they will cancel in the nuclear modification factors, as has been observed in an HEF based approach in [13,14], but nevertheless they could smear the saturation effects, depending on which contributes the most. We plan to tackle those interesting studies in the future. and (4.2): .
By performing a change of variables k t → k t t we get: One partial integration leads to the final result for x 2 G 1 (x 2 , k t ): Similar calculations give the results (4.5)-(4.9).

B High-k t limit of the gluon distributions in the MV model
The large-k t behavior of the dipole and WW distributions has been derived before (see for instance [52] and [53]). We briefly recall these results and then we calculate F We are interested in the region of large transverse momentum, or equivalently in small values of dipole sizes r. Therefore, we expand S F to first non-trivial order in r 2 Q 2 s : 2) The first term in the expansion formally gives a delta function δ (2) (k t ). This term will contribute only for values of k t around zero, and not in the region of large k t that is considered here, so it can be safely dropped. The next terms give: 3) The dipole gluon distribution from Eq. (3.2) is: Similarly, to get the perturbative behavior of the WW density, we start from Eq. (3.4), and we expand the adjoint dipole, S A (x, r) 1− r 2 Q 2 s 2 log 1 rΛ +O r 4 Q 4 s log 2 1 rΛ . For x 2 G (1) we get: Using the above results we can calculate the perturbative expansion of the rest of the distributions. For F (2) qg we have: To this order x 2 G (1) (x 2 , k t ) = x 2 G (2) (x 2 , k t ) and Similarly: and