Numerical analysis of the unintegrated double gluon distribution

We present detailed numerical analysis of the unintegrated double gluon distribution which includes the dependence on the transverse momenta of partons. The unintegrated double gluon distribution was obtained following the Kimber-Martin-Ryskin method as a convolution of the perturbative gluon splitting function with the collinear integrated double gluon distribution and the Sudakov form factors. We analyze the dependence on the transverse momenta, longitudinal momentum fractions and hard scales. We find that the unintegrated gluon distribution factorizes into a product of two single unintegrated gluon distributions in the region of small values of x, provided the splitting contribution is included and the momentum sum rule is satisfied.


Introduction
One of the most intriguing features observed in QCD is the strong increase of the gluon density with the decrease of the parton longitudinal momentum fraction x. This rise, originally discovered at HERA electron(positron)-proton collider [1][2][3][4], manifests itself as the rise of the cross sections for variety of processes in hadronic collisions when the small values of longitudinal fractions of the partons participating in the scattering are probed. The other consequence of this increase in the gluon density is the increase of the probability of multiparton interactions in hadron-hadron collisions. These are events when in one encounter of hadrons, more than one elementary partonic scattering occurs. These type of processes were originally observed by the AFS collaboration at CERN [5], followed by the measurements at the Tevatron [6][7][8] and more recently at the Large Hadron Collider (LHC) [9][10][11].
The theoretical framework for the description of the multiparton interactions in QCD relies on the assumption of the factorization of the hard double parton interaction in the presence of sufficiently large scales, see for example [12,13], and [14,15] for a recent overview. The case of double-parton scattering is usually described within the framework of the standard collinear perturbative QCD by means of the double parton distribution functions (DPDFs) [13,. The DPDF distributions satisfy DGLAP-type equations [16,17,20,21,26,43]. Similar type of equations for double parton correlations were considered earlier [44,45] in the context of the multiparton correlation functions within jets. These equations contain two terms: a homogeneous term describing the independent evolution of the two chains of partons in the double parton distributions and JHEP01(2018)141 a non-homogeneous term which originates from the perturbative splitting of two partons. The latter contribution is driven by the single parton distribution functions. In principle, some part of the splitting contribution can coincide with the higher order contribution to the single parton interaction, and thus one must perform special subtraction and matching in order to avoid potential double counting when computing the double parton scattering cross section with this term [46].
The collinear framework works well for sufficiently inclusive quantities, but can miss important information about the kinematics, see for example [47]. The more detailed information about the kinematics of the process can be incorporated by using the unintegrated parton distribution functions (UPDFs) or transverse momentum dependent distributions (TMDs) [48][49][50]. The unintegrated parton distribution functions in the case of the single scatterings are widely used in the literature, and naturally occur in the high energy or small x limit of QCD, see for example [51]. It is thus of great phenomenological and theoretical interest to develop the formalism which includes the transverse momentum dependence in the context of the multiparton distributions. There has been recently great theoretical effort to develop a consistent formalism for the double parton distribution functions which include transverse momentum dependence, see [52,53] and [54]. The framework presented in [52,53] is based on the extension of the TMD factorization to the double scattering case for colorless final states.
In the present work we perform detailed numerical analysis of the unintegrated double gluon distribution within the formalism which was proposed in our earlier work [54]. In that work, the double parton distributions with the transverse momentum dependence were constructed from the integrated double parton distributions through the convolution with the Altarelli-Parisi splitting functions and by including the Sudakov form factors. This is an extension of the formalism proposed by Kimber, Martin and Ryskin [55][56][57] for the case of the single parton distributions, which has been successfully applied to a wide range of processes which are sensitive to the transverse momentum dependence.
We perform a thorough numerical analysis of the unintegrated double gluon distribution, in particular the dependence on the transverse momenta and longitudinal momentum fractions of gluons, and also on hard scales. In addition, we analyze non-perturbative contributions in which at least one parton is almost collinear with the parent hadron with a low, non-perturbative value of the transverse momentum. We observe that the distribution in the transverse momentum of the double gluon density shifts towards higher values for increasing values of hard scales and decreasing values of longitudinal momentum fraction x. We also investigate the factorization of the unintegrated double gluon distribution into a product of two unintegrated single gluon distributions. We show that factorization holds for sufficiently low values of x provided the splitting contribution in the evolution equations for the integrated double gluon distribution is included and the momentum sum rule for this distribution holds.
The outline of the paper is the following. In section 2 we recall the evolution equations for the integrated double parton distribution functions. In section 3 we discuss the choice of the initial conditions and the issues related to the momentum sum rules. In section 4 we present numerical results on various aspects of the dependence of the unintegrated double JHEP01(2018)141 gluon distribution on kinematic variables while in section 5 we discuss the factorization issues. Finally, in section 6 we present summary and outlook of our analysis.

Evolution equations for the double gluon distribution
We shall start by recalling the evolution equations for the collinear double parton distribution functions. These equations have been first proposed in [44,45] in the context of jet physics for the perturbative QCD description of the jet structure and later derived in [16,17,20] for the initial state double parton distribution functions. We shall review first the evolution equations for the integrated double parton distribution functions, following results of ref. [24]. We note that the DPDFs also depend on the transverse momentum vector, q ⊥ , which we set to zero. This means that in the Fourier space one integrates over the relative position of the two partons b, see [25,53]. In principle this dependence on q ⊥ could be reinstated in the presented formalism through the inclusion of the appropriate form factor, see for example [27]. Precise description of the transverse distribution of partons and the ensuing correlations between the partons in transverse plane is an outstanding and challenging problem [58][59][60][61]. For now, we shall postpone this problem to a further study and consider distributions integrated over b.
For q ⊥ = 0, the DPDFs in the lowest order approximation are probabilities to find two partons with longitudinal momentum fractions x 1 and x 2 , when probed at two different scales Q 1 and Q 2 [13]. The integrated double distributions are denoted from now on by where a 1 and a 2 refer either to the quark flavor or gluon g. The solution to the two scale evolution of the DPDFs can be cast in the following form [24] D a 1 a 2 = a ,a where Q 2 min = min{Q 2 1 , Q 2 2 }, Q 0 is an initial scale for the evolution and the integration limits take into account kinematic constraints x 1 , x 2 > 0 and x 1 + x 2 ≤ 1.
The functions E ab are parton-to-parton evolution distributions which obey the DGLAP evolution equation, with the initial condition E ab (x, µ 0 , µ 0 ) = δ ab δ(1−x). The functions E ab have the interpretation of the Green's functions for the DGLAP evolution, and using them we can construct the evolved single PDF as follows eq. (2.2) is a sum of two contributions which are schematically shown in figure 1 (2.5) The first, homogenous term, D (h) a 1 a 2 , is proportional to the double parton distribution and corresponds to the independent evolution of two partons from the initial scale Q 0 to Q 1 and from Q 0 to Q 2 . The second, non-homogeneous term, D (nh) a 1 a 2 , contains the distribution which describes the splitting of the parton a → a a . Notice on the right hand side of eq. (2.6) the single PDFs, D a , taken at the splitting scale Q s along with the real emission leading order (LO) Altarelli-Parisi splitting functions, In the LO approximation, the flavor of the second parton, a , is uniquely determined from the splitting a → a . The single distributions D a are evaluated at (x 1 + x 2 ) due to conservation of the parton longitudinal momentum in the evolution.

Initial conditions for the evolution equations
Similarly to the case of the DGLAP evolution equations for single PDFs, the evolution equations for the double parton distributions need to be supplemented by the appropriate non-perturbative initial conditions at a scale Q 0 , usually of the order of 1 GeV. The initial conditions for the single PDFs are usually parametrized by some flexible functional

JHEP01(2018)141
forms with many free parameters fixed by fits to the experimental data. There are a few constraints that need to be obeyed by the initial conditions, and these are momentum and quark sum rules. For the double distributions, such sum rules relate the integrals of the DPDFs to the single PDFs, thus placing important constraints on the allowed forms of the initial conditions. This makes the construction of the initial conditions for the DPDFs much more complicated, see [62][63][64] for various proposals. For the case when Q 1 = Q 2 = Q, the momentum sum rules for the DPDFs can be expressed in the following form where D a 1 (x 1 , Q) is the single parton distribution function, known from global fits to hard scattering data. The above momentum sum rule can be interpreted in the following way: the ratio D a 1 a 2 /D a 1 can be viewed as the conditional probability for finding parton a 1 while the longitudinal momentum x 2 of the second parton a 2 is fixed. Therefore the total momentum carried by the remaining partons (except parton a 2 ) is equal to (1 − x 2 ). The quark number sum rule is of a similar form, which can be found in e.g. [22]. In this paper we restrict ourselves to the single channel with gluon only, but the entire analysis could be extended to include quarks as well. In such a case, we are only concerned about the momentum sum rule (3.1) with a 1 = a 2 = g. Assuming that the double gluon distribution is a symmetric function of parton momenta fractions, D gg (x 1 , x 2 ) = D gg (x 2 , x 1 ), it can easily be derived that the rule (3.1) is also valid for the second variable.
In the following, we shall use the initial conditions suggested by the construction proposed in [64] which leads to the initial conditions for both the single and double gluon distribution satisfying the momentum sum rules. The construction is based on the observation that the appropriate set of functions for the initial conditions can be chosen from the set of Dirichlet distributions. Taking advantage of a particular form of the single gluon distribution at the initial scale Q 0 = 1 GeV, given in [65], one can demonstrate that the double gluon distribution obeys the momentum sum rule (3.1), see [64] for details. Note that, the initial double gluon distribution (3.3) is completely determined by the parameters of the single gluon distributions, N k g , α k g , β k g , known from the global fits [65]. It is also worth emphasizing that the distribution (3.3) is not a product of two single gluon distributions, even for small values of parton momentum fractions.
The momentum sum rule is preserved by evolution equations for the double parton distribution [16,17]. Thus, the rule (3.1) is independent of the hard scale Q ≡ Q 1 = Q 2 at JHEP01(2018)141 which the double gluon distribution is defined. The distribution (3.3) will be used in the forthcoming discussion as the initial condition for the evolution equations at some initial scale Q 0 and we shall refer to it as the GBLS 3 input. For the comparison, we also use the distribution proposed in [22] (referred to subsequently as the GS input), This form is approximately factorizable into product of two single PDFs for small values of momentum fractions and includes the correlating factor which guarantees smooth vanishing of the distribution when (x 1 + x 2 ) → 1. This ansatz does not satisfy the momentum sum rule (3.1), although we checked that the violation is relatively small for almost all values of x, that is below 5% percent for x < 0.5, and only becomes significant, up to 30% level, for higher values of x.

Unintegrated double gluon distribution
In this section we present the construction of the double gluon distribution which is unintegrated over the transverse momenta. The framework discussed in the current work, proposed recently in [54], allows to construct the unintegrated double parton distributions which additionally depend on parton transverse momenta, k 1⊥ and k 2⊥ . In the following analysis, we are only interested in the unintegrated double gluon distribution, , though similar investigations can be performed for the other distributions. 1 In the proposed construction [54], which follows the original proposal of Kimber-Martin-Ryskin for the single PDFs [55,56], the dependence on the transverse momenta is obtained from unfolding the last step in the DGLAP-like evolution equations (2.2) reduced to the gluonic sector. The integrated double gluon distribution is the sum of two contributions described in the previous section and illustrated in figure 1, Similarly, the unintegrated double gluon distribution can also be written as a sum of the homogenous and non-homogeneous contributions,

Homogeneous contribution
The form of the homogeneous contribution, f gg , in the region where transverse momenta are perturbative, k 1⊥ , k 2⊥ > Q 0 , is given by [54]

JHEP01(2018)141
In the above equation, T g is the well known Sudakov form factor which sums virtual gluon emissions, where P gg is the real part of Altarelli-Parisi gluon-gluon splitting function in the leading order approximation, The integration limits in eq. (4.3) impose the following condition for the longitudinal momentum fractions which replaces the condition for the integrated double parton distributions, 0 < x 1 +x 2 < 1.
The parameters ∆ in eqs. (4.3) and (4.4) regularize the integrals at z = 1 and could depend on transverse momenta of partons, k 1⊥ and k 2⊥ for real emission and p ⊥ for virtual one, as well as on the hard scales, Q 1 and Q 2 . Thus, eq. (4.6) introduces non-trivial correlations between transverse and longitudinal momenta among both partons. We consider two options in our numerical studies for the kinematics of parton emissions.
1. Strong Ordering (SO): the strong ordering of transverse momenta of emitted partons, typical for DGLAP dynamics, leads to In this case, the transverse momenta are bounded from the above by the hard scale, k i⊥ < Q i .
2. Angular Ordering (AO): the angular ordering of parton emissions (AO) due to the coherence effects in spacelike cascades [66][67][68][69] leads to another cutoff [56,57] In this case the transverse momenta are no longer bounded by the scales Q i .
In figure 2 we show the results for the homogeneous unintegrated double gluon distribution as a function of the transverse momentum k 2 ⊥ ≡ k 2 1⊥ = k 2 2⊥ for four pairs of (x 1 , x 2 ) and fixed hard scales Q 2 = Q 2 1 = Q 2 2 = 100 GeV 2 . We compare the evolution which starts from two discussed initial conditions, given by eqs. (3.3) and (3.4). The calculations were performed for the cutoffs ∆ i which were taken according to the strong ordering (SO) scenario (4.7). We observe that the results for the different inputs are very similar to each other. This is somewhat surprising since the two inputs are rather different. The largest discrepancy is observed at the largest values of x 1 . This is the region where the effects due The distribution , for Q 2 = 100 GeV 2 and the indicated values of (x 1 , x 2 ), starting from the GBLS 3 (3.3) (blue) and GS (3.4) (red) inputs. The cutoffs ∆ i were imposed according to the strong ordering condition (4.7).
to the constraint from the momentum sum rule are most relevant. But even that difference is not large, and this is partially due to the fact that the GS input approximately satisfies the momentum sum rule, up to few percent for the relevant region in x 1 , x 2 . We also see that the unintegrated double gluon distribution vanishes at k ⊥ = Q, as expected in the SO case, and the suppression at the lowest values of k T is due to the Sudakov form factor.
In figure 3 we investigate the effect of different values of Q 2 on the unintegrated double gluon distribution for four fixed values of (x 1 , x 2 ) in the strong ordering case (4.7). We observe that the distribution in k ⊥ is shifted towards higher transverse momenta with increasing values of the hard scale Q 2 , which leads to the greater spread in k ⊥ . Also the peak of the distribution at fixed Q 2 is shifted towards larger values of transverse momenta with decreasing (x 1 , x 2 ). This is due to the fact that as the gluon momentum fraction x decreases, the probability of the gluon splitting increases. The suppression at low values of k ⊥ is increasing with Q 2 due to the Sudakov form factor. The visible sharp cutoff in transverse momenta at the value of the hard scale is consistent with the strong ordering cutoff.
In figure 4 we compare the homogeneous unintegrated gluon distributions in the strong ordering (4.7) and the angular ordering (4.8) cases, plotted as a function of the transverse momentum k 2 ⊥ = k 2 1⊥ = k 2 2⊥ for fixed hard scales Q 2 ≡ Q 2 1 = Q 2 2 = 100 GeV 2 . We should stress at this point that the angular ordering is not implemented in the strict sense since  this condition is only incorporated in the last step of the evolution rather than in the entire cascade, like in the original formulation for single gluon distribution [67][68][69]. Nevertheless, we compare this model with the strong ordering case to illustrate the potential impact of the terms beyond the strong ordering scenario. As expected, we see that the solution with the angular ordering exhibits tails in transverse momentum which extend well beyond the hard scale Q. The tails are most prominent for very small values of x. This can be easily understood since one needs to have x < 1 − ∆ which for angular ordering leads to the condition k ⊥ < Q(1/x − 1) as compared to the condition k ⊥ < Q(1 − x) for the strong ordering. Therefore, at low values of x and high k ⊥ , angular ordering allows for more phase space for emissions. At low values of k ⊥ , the calculation with angular ordering is typically lower than the calculation using strong ordering condition which stems from the fact that the suppression from the Sudakov form factor is larger for the angular ordering condition.

Contributions from non-perturbative regions
In ref. [54] we also provided formulae for the unintegrated double parton distributions in which at least one of the parton transverse momentum is below the initial scale, i.e. for k 1⊥ > Q 0 and k 2⊥ ≤ Q 0 , gg (x 1 , x 2 , k ⊥ , k ⊥ , Q, Q) from eq. (4.3) as a function of the transverse momentum k 2 ⊥ ≡ k 2 1⊥ = k 2 2⊥ in the strong ordering (4.7) (blue curve) and the angular ordering (4.8) (red curve) cases, for the indicated values of (x 1 , x 2 ) and the GBLS 3 input.
Notice that in such a case f (h) gg only depends on perturbative k 1⊥ since k 2⊥ is integrated out in the nonperturbative region below Q 0 . This is reflected in the dependence of D (h) gg on Q 0 . In the opposite case, k 1⊥ ≤ Q 0 and k 2⊥ > Q 0 , we have Finally, for the non-perturbative region, k 1⊥ ≤ Q 0 and k 2⊥ ≤ Q 0 , i.e. there is no transverse momentum dependence, since the transverse momenta have been integrated out below the initial scale Q 0 and the function (4.11) is proportional to the integrated non-perturbative double parton distribution at the initial scales.
In figure 5, we show the non-perturbative gluon contributions given by eq. (4.9), which correspond to the configurations in which one gluon has the perturbative transverse momentum above the cutoff Q 0 , while the other one is non-perturbative. In such a case, the double gluon distribution depends only on transverse momentum of the first gluon. We see that generally the contributions from the non-perturbative regions are numerically much smaller than the perturbative ones, though they can become important for small scales. gg (x 1 , x 2 , k 1⊥ , Q, Q) from eq. (4.9) as a function of k 2 ⊥ = k 2 1⊥ for Q 2 = 10, 100 GeV 2 and the indicated values of (x 1 , x 2 ).

Non-homogeneous (parton splitting) contribution
In ref. [54] we also considered the contribution to the unintegrated double parton distributions originating from a splitting of a single parton into a pair of partons, f (nh) This contribution is generated by the special solution of the non-homogeneous evolution equations for the integrated double parton distributions, D (nh) In the pure gluon case, the integrated double gluon distribution due to this contribution is given by gg under the integral, discussed at length in [54], corresponds to the perturbative gluon splitting g → gg at the scale Q s ,   Figure 8. The homogeneous (4.3) (upper plots) and non-homogeneous (4.15) (lower plots) distributions in the strong ordering case as a function of (k 2 1⊥ , k 2 2⊥ ) for the indicated values of (x 1 , x 2 ) and Q 2 . Contour projections (with 20 equidistant contours) on the transverse momentum plane are shown on the right.
gluon splitting, In the strong ordering case (4.7), the transverse momenta obey the condition k i⊥ < Q i . Thus, for given k 1⊥ and k 2⊥ , the upper limit in the integral over Q 2 s in (4.14) is given by Q 2 min = min{k 2 1⊥ , k 2 2⊥ }. In such a case, the step functions in (4.14) are automatically

JHEP01(2018)141
satisfied, which allows to rewrite the above formula with the help of the distribution (4.12), In the angular ordering case (4.8), the transverse momenta are not bounded and we have to use formula (4.14) together with (4.13) in the numerical analysis. This would mean that when k 1⊥ , k 2⊥ ≥ Q 1 , Q 2 , the splitting could only occur at the scale provided by Q 2 min = min{Q 2 1 , Q 2 2 }, and then the region up to k 1⊥ , k 2⊥ would evolve as two separate parton ladders. However, the parton splitting which generates two partonic ladders can occur at any point of the evolution, and therefore such a situation is unphysical. Since the transverse momentum dependence in the framework considered is anyway generated in the last splitting by convoluting with the splitting functions, therefore we shall argue that also in the case of the angular ordering we should also use Q 2 min = min{k 2 1⊥ , k 2 2⊥ } in eq. (4.14). In figure 6, we show the relative size of the homogeneous (4.3) and nonhomogenous (4.15) contributions by showing them as a function of the transverse momentum k ⊥ ≡ k 1⊥ = k 2⊥ for Q 2 = 10, 100 GeV 2 and the indicated values of (x 1 , x 2 ). We see that the non-homogeneous (splitting) contribution is generally much smaller than the homogenous one. As expected, the significance of the splitting contribution rises with increasing Q 2 .
We also present the unintegrated double gluon distribution for the case with different transverse momenta, k 1⊥ = k 2⊥ . Slices of this distribution for fixed values of k 2⊥ are shown as a function of k 1⊥ in figure 7. The full two-dimensional distributions in (k 1⊥ , k 2⊥ ) are shown in figure 8 together with the contour plots. From the contour plots, it is evident that the distribution originating from the non-homogeneous contribution is more perturbative, i.e. the maximum of the distribution is shifted towards higher transverse momenta and the region of very low transverse momenta is depleted.

Factorization of double gluon distributions
In this subsection we shall investigate to what extent the factorization of the unintegrated double gluon distribution into a product of two single unintegrated gluon distributions, , plotted as a function of k 2 ⊥ for the indicated small values of (x 1 , x 2 ) and Q 2 = 10, 100 GeV 2 , starting from the GBLS 3 input. The breakdown of f gg into the homogeneous (4.3) (dashed blue) and non-homogenous (4.15) (dashed red) components shows that the non-homogeneous component is essential for factorization. eq. (5.1). The results are shown as a function of k ⊥ ≡ k 1⊥ = k 2⊥ for fixed Q 2 and (x 1 , x 2 ). We observe in figure 9 that for small values of (x 1 , x 2 ), the sum of the homogeneous and non-homogeneous contribution is very close to the product of two single gluon distributions. This shows that the factorization of the solution works very well in the small x regime and for large value of hard scale, Q 2 = 100 GeV 2 . It has to be stressed that the factorization works for the sum of the two contributions, gg , but not for the homogeneous term only as might be naively expected. We also observe that the factorization is violated for large values of x 1 or x 2 , see figure 10. This is also seen in figure 11, where we plot the ratio as a function of (x 1 , x 2 ) for fixed transverse momenta, k 2 1⊥ = k 2 2⊥ = 10 GeV 2 and Q 2 = 100 GeV 2 , in the strong ordering case.
The factorization of the distribution f gg is an interesting fact that deserves some explanation. One could naively expect that since the homogeneous term is equivalent to the evolution of two disconnected ladder graphs, the factorization would occur with this term    Figure 11. The ratio (5.2) as a function of the longitudinal momentum fractions (x 1 , x 2 ) for fixed values of transverse momenta, k 2 1⊥ = k 2 2⊥ = 10 GeV 2 , and Q 2 = 100 GeV 2 , in the strong ordering case.
only. On the contrary, it is the sum of the homogenous and non-homogenous contributions which is equal to the product of two single gluon distributions, f To understand this feature, one needs to investigate the factorization of the integrated double gluon distribution, D gg (x 1 , x 2 , Q 1 , Q 2 ), since we expect that the factorization at this level leads to the factorization of the unintegrated double gluon distribution.  Figure 12. The comparison of the integrated gluon distribution x 1 x 2 D gg (x 1 , x 2 , Q, Q) (dashed black line) with the product of two single gluon distributions, x 1 D g (x 1 , Q)x 2 D g (x 2 , Q), (solid magenta) as a function of Q 2 for the indicated values of (x 1 , x 2 ) and the GBLS 3 input. The breakdown of D gg (x 1 , x 2 , Q, Q) into the homogeneous (dashed blue) and non-homogenous (dashed red) components, see eq. (4.1), shows that the non-homogeneous component is essential for the factorization.
In figure 12, we show the integrated double gluon distribution D gg as a function of the hard scale Q = Q 1 = Q 2 for fixed values of (x 1 , x 2 ). We indeed observe that the sum of the homogenous and non-homogenous terms nicely factorizes into a product of two single integrated gluon distributions, provided x 1 and x 2 are sufficiently small and Q 2 is large. The latter condition is related to the fact that at low Q 2 , the solution is still close to the form of the initial condition which does not factorize in the case of the GBLS 3 input. So even though the initial condition does not factorize at small x, after the evolution to large scales the solution becomes factorized for small x.
In order to further investigate the origin of the factorization, we also performed calculations for the GS initial conditions. We found that in this case the factorization at large scales holds to a lesser degree than for the GBLS 3 input, though approximately it is still valid. Finally, we also tested initial conditions which violate the momentum sum rule more strongly than the GS input [22], and we found that in this case the factorization of the resulting solution is stronger violated.
The crucial difference between the analyzed inputs is the fact whether or not they satisfy the momentum sum rule. The best quality factorization at small x occurs for the GBLS 3 input which satisfies the momentum sum rule exactly by construction. On the other JHEP01(2018)141 hand, the GS input does not satisfy the sum rule, but the violation is of the order of a few percent (for x < 0.5), which translates to a similar in size effect for the factorization. This suggests that whether or not the factorization at small x and large Q 2 holds is related to the fulfillment of the momentum sum rule for the double PDFs. We stress, that this is true even if there is no factorization in the initial condition. However, more detailed analytic studies are necessary to quantify this observation, which we postpone for the future.

Summary and outlook
We performed detailed numerical analysis of the unintegrated double gluon distribution with the transverse momentum dependence. The construction of such distributions, presented in [54], is based on the extension of a similar formalism for the construction of the single unintegrated parton distributions [55,56]. The unintegrated double parton distributions in [54] were obtained as a convolution of the integrated parton distributions with the Altarelli-Parisi splitting functions and the Sudakov form factors.
As a general conclusion from the presented analysis, the transverse momentum dependent double gluon distribution is shifted towards larger values of transverse momenta with increasing values of the hard scales Q 1,2 and decreasing values of the longitudinal momentum fractions x 1,2 . The model with angular ordering in the last step of the evolution leads to a distribution which extends further in transverse momenta than in the strong ordering scenario in which the transverse momenta are cut off by the hard scales. The homogeneous contribution dominates over the non-homogeneous one, for all the values of transverse momenta, at least for the type of initial conditions considered in this work. We observe that the non-homogeneous contribution, which originates from the splitting, is more perturbative being shifted towards higher values of transverse momenta than the homogeneous contribution. The sum of the homogeneous and non-homogeneous contribution factorizes into the product of two single unintegrated gluon distributions at small values of x 1,2 and high hard scales. We showed that this factorization stems from the factorization of the underlying integrated parton distributions and is contingent upon the validity of the momentum sum rule for the initial conditions.
There are number of possible avenues that could be further explored. First of all, we have only considered the distributions for q ⊥ = 0, i.e. integrated over the relative transverse positions of the two gluons. Full dependence on q ⊥ should be included when one wants to use these distributions for the evaluation of the double parton scattering cross sections. In the simplest scenario, this could be done by including the appropriate form factors. However, it is possible that there will be large correlations between the transverse momenta of the partons and the momentum transfer q ⊥ , see [15,46,53] for recent studies. Second, the presented analysis was performed for the gluon sector only. The extension to include quarks can be done, and in principle does not present any new technical difficulties. The biggest challenge is the appropriate choice of the initial conditions for the evolution which would satisfy both the momentum and the quark number sum rules simultaneously. So far this problem has not been fully solved. Finally, it would be also interesting to further explore the relation between the momentum sum rule and the factorization of the double distributions for small values of the longitudinal momentum fractions.