Drell–Yan production with the CCFM-K evolution

We discuss the Drell–Yan dilepton production using the transverse momentum dependent parton distributions evolved with the Catani–Ciafaloni–Fiorani–Marchesini–Kwieciński (CCFM-K) equations in the single loop approximation. Such equations are obtained assuming angular ordering of emitted partons (coherence) for x∼1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x\sim 1$$\end{document} and transverse momentum ordering for x≪1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x \ll 1$$\end{document}. This evolution scheme also contains the Collins–Soper–Sterman (CSS) soft gluon resummation. We make a comparison with a broad class of data on transverse momentum spectra of low mass Drell–Yan dileptons.


Introduction
The Drell-Yan (DY) dilepton production [1] is one of the most intensively studied processes in particle physics. The existence of a hard electroweak probe (photon or Z boson), which does not interact strongly and decays into a pair of leptons, provides a clear experimental signature of the partonic interactions in the colliding hadrons and significantly simplifies their theoretical description. For these reasons, the DY process is very efficient tool for the investigation of hadronic structure [2], in particular the distributions of partons' transverse momentum. The key concept of the theoretical analysis of DY scattering, called factorization, is a separation between long-range and short-range degrees of freedom. The basic, collinear factorization theorem assumes no transverse momenta of partons in hadrons [3]. Its application to the Drell-Yan production is very well established and commonly used. The present state of art calculation for the DY process includes next-to-next-to-leading order (NNLO) QCD corrections [4][5][6].
There are several kinematical regimes where the fixedorder collinear QCD does not provide good description of data. In particular, when transverse momentum of lepton pair is much smaller than it's invariant mass, q T M, the large logarithms log n (M/q T ) occur in all orders of perturbative expansion. These corrections are effectively resumed within the Collins, Soper and Sterman (CSS) approach [7]. As a result, the collinear factorization need to be replaced by a new factorization based on the transverse momentum distributions (TMDs) [8]. For state of art TMD analyses of Drell-Yan process, see [9][10][11][12][13][14]. It is interesting to ask what is the transition from the small transverse momentum region to the region of q T ∼ M, where fixed order perturbative QCD should apply. As it was shown in [15], when one considers moderate values of M ∼ 5 − 19 GeV, the fixed-order predictions underestimate data in the region of q T ∼ M.
In this paper, we examine in detail the low mass DY productions using the k T -factorization approach combined with the transverse momentum parton distributions defined in the Catani-Ciafaloni Fiorani-Marchesini branching scheme (CCFM) [47][48][49][50]. The original motivation for the CCFM branching was to extend angular ordering (coherence) of soft parton emission in the space-like branching from the region of x ∼ 1 to the region of small x 1. In this way, a unified evolution equation for transverse momentum dependent gluon distribution was found with angular ordering in both regions of x in the approximation called all loop. In the fully inclusive case, this equation interpolates between the DGLAP equation at moderate x and the BFKL equation obtained in the small x limit. It is worth emphasizing that the CCFM branching scheme for x ∼ 1 contains the CSS resummation of soft gluon emissions [51]. The Monte Carlo implementation of the all loop CCFM branching scheme was done in [52][53][54]. A general Monte Carlo scheme for QCD evolution was also constructed with the Parton Branching method [55][56][57] and subsequent analyses were presented in [58,59].
In the region of large or moderate values of x, when the small x coherence can be neglected, the CCFM scheme gives the DGLAP evolution equation with coherence at large x only. This approximation, called single loop, was studied in [50,52] for the gluon distribution. The extension of this scheme in terms of evolution equations for both quark and gluon distributions was proposed by Kwieciński in [60]. This is why we call them the CCFM-K equations. The parton distributions which are obtained by solving these equations depend on transverse momenta and their properties were analyzed in [61][62][63][64][65]. The first analysis of the weak gauge boson production with the CCFM-K equations was done in [66] while the low mass DY production with photon was addressed in [67]. Similarly to the all-loop case, the one-loop CCFM-K equations contain a part of the CSS resummation of soft parton emissions [63].
The main goal of the presented analysis is a comprehensive analysis of all available data on transverse momentum spectra in low mass DY production with the CCFM-K evolved parton distributions and the leading order cross sections computed in k T factorization. We assume the most eco-nomical form of the initial conditions for the evolution with only one adjustable parameter. In this way, we concentrate on the most important effects of the CCFM-K evolution which are responsible for good description of data for q T ∼ M. The small q T description is acceptable in most cases, especially for higher masses, although precise comparison should involve more adjustable parameters in the initial conditions like in the CSS approach. This is left for future studies.
The paper is organized as follows. In Sect. 2 we give an overview of the CCFM framework and its version proposed by Kwieciński in which quark and gluon transverse momentum dependent distributions and the CCFM-K evolution equations are introduced. We also discuss the relation of the CCFM-K approach to the CSS formalism (with the full derivation presented in "Appendix A"). In Sect. 3 we describe the application of the discussed formalisms to the leading order DY cross section with both on-shell and offshell matrix elements in the CCFM-K case. In Sect. 4 we show numerical results and compare them with the low mass DY data. We summarize in Sect. 5.

Branching kinematics
Below we describe kinematics of the CCFM parton branching schemes in two approximations -single and all loop. We work in the Sudakov base with two light cone vectors in which the momenta in the branching process shown in Fig. 1 are given by Notice that the momentum fraction are proportional to the plus/minus components, e.g.
The emitted parton momentum can be found from the momentum conservation where x i−1 > x i . Denoting the transverse component by p i T = k (i−1)T −k i T = (0, p i T , 0) and assuming that p 2 i = 0, one can compute the minus component to find The rapidity of the emitted parton is given by where z i = x i /x i−1 < 1. In the massless case y i = − ln tan(θ i /2) where θ i is the emission angle with respect to the z axis defined by the collinear momenta P 1 and P 2 , therefore The CCFM branching scheme [50] is defined with the help of the rescaled transverse momentum of the emitted parton, Thus, the transverse momentum conservation at the vertex i reads From (7) we obtain for the modulus The single loop approximation is defined by the condition thus for z i → 0 we find the transverse momentum ordering of emitted partons For finite z ∼ 1, however, from (10) applied to (11) we have and for z i−1 = (x i−1 /x i−2 ) → 1 we obtain angular ordering of parton emissions Such a phenomenon is called coherence [47][48][49][50]. Thus, in the single loop approximation partons are emitted with transverse momentum ordering for z → 0 and angular ordering for z → 1. The all loop approximation is defined by the condition which is equivalent to the condition giving the angular ordering condition (14). Thus, in the all loop approximation partons are emitted with angular ordering for any value of z.

CCFM equation in all loop approximation
The CCFM branching schemes allow to define the corresponding parton distributions. In the all loop approximation [50] only the gluon distribution f g is defined up till now through the equation which relates the gluon distribution at vertex i with the gluon distribution at vertex (i − 1). In the above, k T = | k T | and q = | q| are transverse momenta depicted in Fig. 1 and S is the Sudakov form factor given by where P gg is the gluon-gluon splitting function (27), which resums virtual corrections for z → 1. For z → 0, the virtual emissions are resumed by the non-Sudakov form factor Notice that only the 1/z part of P gg (z) is present under the integral. The first theta function in Eq. (17) reflects ordering (15) in which Q is a hard scale which terminates the CCFM evolution, Therefore, for given x and √ S, the hard scale determines the maximal emission angle The second theta function in (17) imposes the condition which assures that α s (q) 1 and the CCFM evolution scheme is perturbative. From the third theta function in (17), we find the following condition for the real gluon emission which allows to avoid singularity of P gg (z) at z = 1. Non-perturbative effects are encoded in the initial condition, Equation (17) can be used for the Monte Carlo generation of a parton cascade with angularly ordered emissions which leads to the gluon distribution f g . Intensive studies with the CCFM-K equations in all-loop approximation were done using Monte Carlo program CASCADE [52][53][54].

CCFM-K evolution equations
The mixing between the transverse and longitudinal variables in the theta function θ(Q−zq) prevents writing Eq. (17) in the form of an evolution equation. However, this can be done in the single loop approximation in which the branching scheme leads to the CCFM-Kwieciński (CCFM-K) evolution equations for both quark and gluon distributions. The evolution scale is defined in such a case by the rescaled momentum Q.
In order to obtain the CCFM-K evolution equations in the single loop approximation, the branching conditions in Eq. (17) are replaced by [60][61][62][63] Thus, for z → 0, the angular ordering is replaced by the transverse momentum ordering while for z → 1 the angular ordering is still valid. In addition, quark splittings, q → qg, q →qg and g → qq, are taken into account which allow to introduce quark distributions f i=1,...,2N f in addition to the gluon distribution f g . In this way we obtain [63] f i (x, k T , Q) where the argument of the parton distributions on the r.h.s. equals The one loop real emission splitting functions are given by 1 where C F = 4/3, C A = 3, T R = 1/2 and N f is the number of active flavours. Notice that after integrating both sides of Eq. (25) over k T , the ordinary DGLAP equations are found for the collinear quark and gluon distributions, Equation (25) can be written with the help of the Fourier transformatioñ which for b = 0 gives the PDFs (28). Since the parton distributions depend on k T = | k T |, we perform the azimuthal angle integration with help of the relation and obtain the parton distributions which depend on b = | b|, Thus, taking the Fourier transform of both sides of Eq. (25), we find the evolution equations which are diagonal in b: These are the CCFM-K equations which we use in our forthcoming analysis. As expected, for b = 0 we obtain the DGLAP evolution equations for the collinear PDFs (28), i.e.
It should be emphasized that the studies with the CCFM-K equations were also done using the Parton Branching (PB) method for the construction of the TMD parton distributions [55,56] which is based on Monte Carlo algorithms. Recently, the low mass DY production was analyzed with this method in [59]. The main difference between our approach and the PB method, aside from technical issues, lies in the treatment of the strong coupling constant α s in the CCFM-K equations. We keep it outside the integrals on the rhs of Eq. (31) with the scale given by the evolution variable Q, whereas in the PB method α s is inside the integrals over z since it depends on the transverse momentum of an emitted parton, In such a case, a cutoff on the upper limit of z is necessary to avoid the Landau pole in α s (k T ).

Initial conditions and b-dependence
In order to solve Eq. (32), we need initial conditions specified as functions of x and b at some perturbative scale Q 0 QC D . They have to fulfill the conditions saying that for b = 0 the collinear PDFs are recovered. Thus the simplest possible choice is given in the factorized form where q i and g are the LO collinear quark and gluon distributions at scale Q 0 and the non-perturbative form factor obeys the condition F(0) = 1. In the forthcoming analysis we will use the gaussian form factor with one free parameter b 0 , In principle, different form factors can be used for quarks and gluons. However, with the common form factor, it is possible to write the solution of the CCFM-K equations for any value of Q 2 as a product where f CCFMK i,g is the solution for F(b) ≡ 1. This is because the Eq. (32) are homogeneous, thus the multiplication by the common form factor F(b) can be done at the beginning or the end of the evolution. In this way, the perturbative and non-perturbative dependences of the solution are clearly separated.
This effect is shown in Fig. 2 for the singlet quark distribution, f S = i f i , plotted as a function of b for fixed x = 10 −2 . The dashed curves are the initial conditions at Q 0 = 1 GeV with the MSTW08 LO PDFs [68] while the solid curves are evolved to Q = 100 GeV. On the left plot, the b-dependence of the evolved curve is purely perturbative while on the right plot the curves were multiplied by the form factor (35). We see that its impact is the strongest for large values of b while for small values, the b-dependence of the full solution remains perturbative. After the Fourier transformation to the k T -space, we find broadening of the parton distributions due to the CCFM-K evolution, studied in detail in [62,64,65].

Relation to CSS resummation
The CCFM-K equations contain a part of the Collins, Soper and Sterman (CSS) resummation of soft parton emissions in the limit z → 1. In "Appendix A", we present the proof that for the values of the parameter b such that the solution to the CCFM-K equations (32) is given by the CSS formulas [7] where the parameters A (1) q,g and B (1) q,g are defined in Eq. (90). The above formulas were derived by picking large logarithms, ln(Qb) and ln(1/Q 0 b), in the CCFM solution. It should be noted, however, that Eq. (38) do not contain the NLL (next-to-leading logarithmic) terms proportional to α 2 s under the integrals (see Sect. 3.3) since the splitting functions in Eq. (32) are in the leading order approximation.
In Fig. 3 we present the numerical solution to the CCFM-K equations (32) with F(b) = 1 (solid lines) against the CSS approximation (38) (dot-dashed lines) for the singlet quark (left plot) and gluon (right plot) distributions. For the chosen scales, Q 0 = 1 GeV and Q = 100 GeV, the range (37) corresponds to b ∈ [10 −2 , 1] GeV −1 . We see that the CSS approximations extracted from the CCFM-K equation works reasonable well for b ∈ [10 −2 , 10 −1 ] GeV −1 . For b = 10 −2 GeV −1 the two analyzed curves coincide, which results from the observation that for the scale Q = 1/b = 100 GeV, corresponding to this point, both the CSS formulas (38) and the CCFM-K solution are equal to the collinear PDFs at the scale Q. This is obvious for Eq. (38), while for the CCFM-K solution it is a manifestation of the DGLAP limit (33) at b = 0, which becomes already effective for b = 10 −2 GeV −1 . Beyond the lower limit in (37), the CSS approximation significantly deteriorates and the approximation (38) cannot describe the CCFM-K solutions.
The condition 1/Q < b which was necessary for us to find the connection between the CCFM-K and CSS approaches is not present in the original CSS formulation [7], where b can be arbitrary small. Nevertheless, recent studies [69] introduces such a constraint, i.e. b is limited from below by b min ∼ 1/Q. Analyzing this idea in the context of the CCFM-K approach, however, would go beyond the main thrust of our analysis.

Drell-Yan cross section with k T -dependent PDFs
The Drell-Yan cross section differential in photon's momentum is given by where (y γ , M, q T ) are photon's rapidity, virtuality and transverse momentum while W μ μ is the trace of the hadronic tensor W μν . With the lowest order matrix element for the process qq → γ * , the trace is given by the transverse momentum factorization formula where k 1T , k 2T are quark transverse momenta, x 1,2 are their longitudinal momenta and f i , f¯i are transverse momentum dependent quark/antiquark distributions taken at the scale Q = M.

On-shell matrix element
The trace in (40) is the squared matrix element of the process q(k 1 )q(k 2 ) → γ * in the lowest order. For the on-shell matrix element, we use the quark/antiquark momenta in the collinear approximation what assures gauge invariance of the matrix elements. In such a case (42) and the DY cross section (39) is given by It is easy to check that after integrating (43) over q T , we find the leading order form of the Drell-Yan cross section with collinear PDFs given by Eq. (28) where σ 0 = 4πα 2 em /3N c . Inserting the delta function to Eq. (43), we find the DY cross section with the b-dependent parton distributions (29) dσ DY on−shell For the parton distributions which depend on b = | b|, the angular integration with the help of relation (30) gives We will use this expression for the comparison with the DY data using the parton distributions which are solutions of the CCFM-K equations with the momentum fractions in the onshell form

Off-shell matrix elements
In approach with the off-shell matrix, the trace (42) is replaced by where μ is the Fadin-Sherman photon-quark vertex [70,71] and the quark/antiquark momenta k 1,2 take into account transverse components They are given by k T i = (0, k T i , 0) for i = 1, 2 while the momentum fractions x i are determined from the momentum conservation at the vertex, q 2 = (k 1 + k 2 ) 2 , which gives It is easy to check that the Fadin-Sherman vertex obeys the gauge invariance relation Computing the trace (49), we obtain where we used momentum conservation at the photon vertex to write the last equality. Notice that because of the transverse mass M ⊥ in the denominator, the off-shell kinematics takes into account the corrections in powers of q 2 ⊥ /M 2 to all orders. With thes results, the cross section (39) is given by Inserting the delta function (45) and performing the Fourier transformation, we obtain the following cross section with the parton distributions which depend on b = | b| where b is the radial part of the two-dimensional Laplacian b =  (47), we see that (56) has different mass dependence, It should be emphasized that the corrections q 2 T /M 2 which are resummed to the factor 1/M 2 ⊥ in (56) are entirely due to off-shellness of the matrix element. In the CSS approach such corrections, if large, signal the breaking of the CSS approximation. However, in the approach with transverse momentum dependent parton distributions (like the CCFM-K approach), they are naturally incorporated in the PDFs and off-shell matrix element. This is the main advantage of this method.
Numerical studies show that the contribution from the terms in the third and fourth lines in (56) is negligible. Therefore, for the same values of x 1 and x 2 , the cross section dσ DY off−shell is suppressed by a factor M 2 /M 2 ⊥ in comparison to dσ DY on−shell . In addition, in the off-shell case the PDFs are taken at larger values of x 1 and x 2 , compare (48) and (52), which additionally suppresses dσ DY off−shell at large q ⊥ . This effect is clearly visible in Fig. 4 where we plot the CCFM-K predictions against the Fermilab R209 data [42]. The solid line corresponds to the on-shell cross section (47) with the CCFM-K parton distributions evolved from the initial conditions (34) at Q 0 = 1 GeV with the MSTW08 LO PDFs [68] and the form factor (35) with b 0 = 2.7 GeV −1 . The dashdotted line is obtained from the off-shell cross section (56) with the same parton distributions. We also show, for general orientation, the prediction from the CSS formula (59) (red dashed line) which is discussed in detail in the next section.

CSS approach
The CSS approach to the DY process has a long history which starts with the pioneering work [7]. In this approach, collinearly colliding quarks emit gluons with a total transverse momentum q T which is balanced by the transverse momentum of the DY boson. The soft and collinear divergences for q T → 0 in real emission are not fully cancelled by virtual corrections and manifest themselves by the presence of large logarithms, log(M/q T ), which are resummed in the CSS approach. This leads for example to evolution equations with two scales in the NNLL approximation. The recent state-of-art analyses of the DY data with the CSS approach up to the order N 3 LL were done in [13,14].
Such a precision in the CSS approach is beyond the scope of our present analysis since we do not aim at a comprehensive description of the data using this approach. It only serves for the comparison with the results of the CCFM-K approach in which q T of the DY boson is the sum of intrinsic transverse momenta of colliding partons, see the formulae (39) and (40). For this reason, we also do not consider the so-called "Y term", which was proposed in [7] to match the CSS formula with the fixed-order results. Thus, we will use the CSS formulae in the NLL approximation with one scale evolution. Nevertheless, all the problems which are encountered in the description of the DY data in this approximation are still present in the analyses done with higher order approximations.
The DY cross section in the CSS approach up to the NLL order is given by where Q = M and x 1,2 = Me ±y γ / √ S. The parton luminosities W iī in the above read where f i/ī are effective quark/antiquark distributions with the MS NLO collinear PDFs q i/ī (z) and g(z) and the coefficient functions The scale μ in (61) is given by μ = c/b * with c = 2e −γ E ≈ 1.12 and TIn this way, b * interpolates between b * = 0 and b * = b max for b → ∞ such that the scale Thus, choosing b max = c/Q 0 , where Q 0 is an initial scale for the DGLAP evolution, we ensure that the collinear PDFs are always defined during the integration over b in (59). In our presentation, we use the MSTW08 NLO PDFs [68] and choose Q 0 = 1 GeV. The power S in the exponent in (60) is given by where the coefficients A q and B q are defined by the general perturbative expansion Introducing B = ln(Q 2 b 2 /c 2 ) and L = L(Q) = ln(Q 2 / 2 ), the LL approximation is defined by the terms proportional to B(B/L) n while in the NLL approximation terms proportional to (B/L) n are added. Thus, in he NLL approximation which we consider, the coefficients are given by [72][73][74] and K = C A ( 67 18 − 1 6 π 2 ) − 10 9 T R N f . By the comparison of the power S given by (65) with that in (38), we see that the CCFM-K equations only partially resum the next-to-leading logarithms since the term proportional to A (2) q , which is formally of the NLL accuracy, is missing in (38). It can be obtained, however, from the CCFM-K equations with the higher order splitting functions. Using the two-loop running coupling constant where β 0 = (11C A − 4T R N f )/12π and β 1 = (153 − 19N f )/24π 2 , and performing the integration in (65), one obtains the final NLL form of S [74], which we use in our presentation The factor W N P = exp{−S N P } in (60) describes the nonperturbative contribution [75][76][77]. In our presentation, we use the form factor from the BLNY fit [78] to the DY data: where a 1 = 0.21 GeV 2 , a 2 = 0.68 GeV 2 ,

Comparison to data
For the comparison with the low mass DY data, we use the CCFM-K approach with both on-shell and off-shell matrix elements (see Sects. 3.1 and 3.2). In this approach, we only have one free parameter, b 0 in the non-perturbative form factor (35), which we somewhat optimized to the value b 0 = 2.7 GeV −1 .
For the initial PDFs we use the MSTW08 LO PDF set [68]. We also show the CSS results at the NLL accuracy with the BLNY form factor (70) and the MSTW08 NLO PDF set [68] (see the previous section). We use this more refined form factor and PDF sets (compared to those of CCFM-K) in order to reach better description of data within the CSS formalism at the NLL accuracy. The results depend to some extent on the value of the parameter b max = c/Q 0 in Eq. (63) but not such that the general conclusions concerning the CSS description should be changed. For example, using Q 0 = 2 GeV makes the curves stronger suppressed for large q T .
Any attempt to have exactly the same set of parameters for both the CSS and CCFM-K approaches leads to significant deterioration of the agreement with the data in one or the other description. This is not a surprise since the CSS and CCFM-K approaches have different starting points (collinear versus k Tfactorization) and are derived in different approximations (NLL versus LL). Therefore, they have to be optimized with respect to the DY data description separately. To check the impact of the choices we made, we produced the results for NLL CSS with the Gaussian form factor (35) (with properly chosen b 0 ) and MSTW08 LO PDF. It turns out that in such a case, the description of data at small q T is worse than for the one with the NLO PDFs and the non-perturbative contribution (70). Moreover, as expected, the rapid fall of the CSS curves at high q T is still present (see also [15] for detailed discussion of difficulties with matching the CSS approach to fixed order calculation and description of data at q T ∼ M).

DY from fixed target experiments
We start with the data from the fixed target experiments E288 [38], E605 [39], E866 [41] and E722 [40]. The cross section Ed 3 σ/d 3 q measured in these experiments is related to (47) and (59) as follows where M is the bin size of the DY pair mass distribution. In addition, the E605, E866 and E722 experiments also measured the cross section in bins of the Feynman variable where the r.h.s. gives the relation between x F and the DY photon rapidity y γ . The energies of the proton projectile were equal to: 200, 300 and 400 GeV (at E288), and 800 GeV (at E605, E866 and E772). These translate into the center of mass energies √ S = 19.4, 23.8, 27.4 and 38.8 GeV, respectively. The experiments differ by the targets used: E288 used Cu or Pt, E605 used Cu, E866 used H or D and E772 used 2 H. All the cross sections were normalized by the number of nucleons in the target nucleus. In what follows, we neglect nuclear effects of the targets and compare unmodified CCFM-K and CCS approaches with such data.
In Fig. 5 we show the data from the E288 experiment [38] for three values of energies and rapidities. At each panel, the transverse momentum dependence of the DY cross section is shown with fixed mass M and rapidity y (or x F ) 2 . The mass range equals M = 4.5−12.5 GeV. In Fig.6 we show the data from the E605 [39] and E866 experiments [41] for √ S = 38.8 GeV, x F = 0.1 and the mass range M = 4.7−15.5 GeV. In Fig. 7 we show the data from the E866 experiment for √ S = 38.8 GeV, three values of x F and the mass range M = 4.7−14.85 GeV. Finally, in Fig. 8 we show the data from the E772 experiment [40] for √ S = 38.8 GeV, 0.1 < x F < 0.3 and the mass range M = 5.5−14.5 GeV.
The data in Figs. 5, 6, 7 and 8 are compared to theoretical curves: CCFM-K on-shell (blue solid curves), CCFM-K offshell (blue dashed-dotted curves) and CSS (red dashed curves). Comparing CCFM-K to CSS we see that at small q T the CSS resummation predicts higher cross-section than CCFM-K and better agrees with the E288 and E605 data. This comes from the fact that the parameters of the non-perturbative form factor (70) were fitted in [78] to the data while in the CCFM-K approach Fig. 5 Transverse momentum dependence of DY cross section: data from fixed target experiment E288 [38] are compared with on-shell CCFM-K (blue solid), off-shell CCFM-K (blue dash-dotted) and CSS (red dashed) approaches Fig. 6 Transverse momentum dependence of DY cross-section: data from fixed target experiments E605 [39] (upper panels) and E866 [41] (lower panels) are compared with on-shell CCFM-K (blue solid), off-shell CCFM-K (blue dash-dotted) and CSS (red dashed) approaches we fixed just one free parameter b 0 in the form factor (35) The motivation for that was to show the potential of the CCFM-K approach to describe the large q T data without going into details of fitting the parameter b 0 , as in the CSS approach. Thus, at larger q T (2-3 GeV, depending on M and x F ), the CSS curves drop rapidly as we do not match them to the fixed order calculation by adding the "Y term". On the other hand, the CCFM-K curves describe the data reasonably well. Note also that for M ∼ 9 GeV, the data from E288 are significantly above the theoretical predictions which is related to the production of ϒ meson, not considered in our calculations.
The E866 and E772 data seems to be systematically above theoretical predictions at small q T , except for a few values of M Fig. 7 Transverse momentum dependence of DY cross-section: data from fixed target experiment E866 [41] are compared with on-shell CCFM-K (blue solid), off-shell CCFM-K (blue dash-dotted) and CSS (red dashed) approaches Fig. 8 Transverse momentum dependence of DY cross-section: data from fixed target experiment E772 [40] are compared with on-shell CCFM-K (blue solid), off-shell CCFM-K (blue dash-dotted) and CSS (red dashed) approaches Fig. 9 Transverse momentum dependence of DY cross-section in proton-proton collisions: data from PHENIX [46] compared with onshell CCFM-K (blue solid), off-shell CCFM-K (blue dash-dotted) and CSS (red dashed) approaches and x F . As before, CCFM-K provides good description of the data at large q T . In general, one sees better description of data for higher DY masses.
Comparing the CCFM-K on-shell and off-shell approaches one sees that the former approach agrees better with data as providing slower drop with q T . The difference is larger for q T ∼ M, as one should expect, see discussion at the end of Sect. 3.2.

DY in proton-proton collisions
We also consider the data from two experiments measuring the DY production in proton-proton collisions at moderate energies: R209 [42] with √ S = 62 GeV and PHENIX [46] with √ S = 200 GeV. For R209 we apply a change of variables, whereas PHENIX is using the cross section Ed 3 σ/d 3 q given by (72). The theoretical results were compared with the R209 data in Fig. 4. CCFM-K provides very good description of data up to q T ∼ 4 GeV and slightly underestimate cross-section for higher q T while CSS overestimates the cross section at small q T and decreases rapidly at high values. For the PHENIX data shown in Fig. 9, CCFM-K gives a better description than CSS, which overestimates the data at moderate values of q T . We note that as for fixed target experiments, the on-shell CCFM-K better describes the data then the off-shell approach.

Summary
Using the CCFM-K parton distributions and the partonic crosssection with on-shell and off-shell matrix element, we analyze the transverse momentum spectra of the DY dileptons from all available low mass data. The overall description of these data is quite good, given the simplicity of the non-perturbative Gaussian input (35) with only one free parameter b 0 , being the Gaussian width. We have chosen optimized value of this parameter for all experiments to show the potential of the CCFM-K approach in the description of the data for large dilepton transverse momentum q T .
However, for small q T we found less successful description, especially in low mass bins. This calls for an approach with more non-perturbative parameters in initial conditions for the CCFM-K evolution akin to the BLNY fit [78] of the non-perturbative form factor (70) in the CSS approach. This is justified since the CCFM-K evolution includes elements of the CSS resummation. In this sense our paper should be treated as a step towards unified description of the low mass DY data, where the CSS approach matched to the fixed order calculation experiences some troubles [15].
One should also note that the presented analysis is based on the LO matrix elements and the CCFM-K evolution equations in the single loop approximation. For these reasons, we decided to postpone the analysis with more complicated non-perturbative input to future studies with the NLO matrix elements and evolution equations. The first attempt in this direction was done recently in [59] using the Parton Branching method. Finally, it is important to stress that the future analysis should also include the Tevatron and LHC data on the weak bosons production.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: The research described in the article is of purely theoretical nature and no experimental data have been collected in the process of working on the paper.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/. Funded by SCOAP 3 .

A Relation of CCFM-K to CSS
Following the method presented in [63], we show the CCFM-K resummation contains the soft gluon resummation of Collins, Soper and Sterman (CSS) [7]. In order to simplify the notation, we apply the Mellin transform to both sides of Eq. (32) to obtain wheref S = 2N f i=1f i is the singlet quark distribution and the Mellin moments of the splitting functions read Notice that for b = 0 in (76) we obtain the DGLAP evolution equations. This fact motivates the following decomposition of the diagonal splitting functions wherē Thus, for b = 0,P qq2 =P gg2 = 0, andP qq1 andP gg1 become the ordinary Altarelli-Parisi splitting functions. This is why we write (76) in the form where for simplicity we suppressed the arguments. We look for the solutions in the form where Inserting (84) to (83), we find We are interested in the approximate solution when b obeys relation (37), i.e. for In such a case, the powers of large logarithms, ln(Qb) and ln(1/Q 0 b), organize the calculations. We will show that in such an approximation the solution is given by the Mellin moments of the CSS formulas (38) where the coefficients Proof The dominant contribution to the integrals with the Bessel function J 0 (u) comes from the region u < 1, therefore, we use the following approximations where θ is the Heaviside step function and c ∼ 1. The precise value of this parameter is important for numerical studies but for simplicity of this analysis we set c = 1. Thus, the quark exponent S q in (85) withP qq2 given by (80) is given by For qb < 1, the argument of the theta function is negative and S q = 0. Thus, we have From the first theta function q > 1/b Q 0 which sets the lower integration limit to 1/b. In this way, we avoid resummation of large logarithms ln(1/Q 0 b) which are shifted to the functions f i,g in (84) and need to be resummed separately. Since 1/b < Q in our approximation, the upper integration limit Q 2 in (93) is not affected by the first theta function. Writing P qq in the form we find the result which agrees with the exponent in (88) where in the last line we neglected terms with subleading powers of logs ln(Qb) after the integration. A similar calculation for S g in (85) leads to the exponent (89) The large logs ln(1/Q 0 b) are resummed using the DGLAP evolution equations. To show this, we write (86) in the integral formf The first integral on the r.h.s. is given by where we used (79). Approximating we notice that in order to get the leading logs ln(1/Q 0 b) we have to assume that μb < 1. In such a case, the integration over z in (98) is not constrained but the integration over μ is limited to μ < 1/b Q. In this way Applying the same approximations to second integral in (97), we obtain From (95), in the integration region, μ < 1/b, we have S q (b, μ) = S g (b, μ) = 0. Therefore, we can set the exponent in I 2 equal one and (97) reads These are the DGLAP equation for the moments of quark the distributions evolved to the scale Q = 1/b. The analogous considerations for the gluon distributions lead to the gluon counterpart of the DGLAP equations. It can easily be checked that for μ > 1/b, the theta function (99) imposes the condition z > 1−1/(μb) leads to subleading logarithms which we neglect.
To make a connection with the collinear PDFs, we note that for the values of b which we consider, b 1/Q 0 ∼ 1 GeV −1 , to good approximation b ≈ 0 in the functionsf i,g in (102). Thus, choosing the initial conditions equal to the Mellin moments of the collinear PDFs, f i (n, b ≈ 0, Q 0 ) = q i (n, Q 0 ), we find the collinear PDFs at the scale 1/b f i (n, b ≈ 0, 1/b) = q i (n, 1/b), f g (n, b ≈ 0, 1/b) =ḡ(n, 1/b).
This concludes the proof that (88) and (89) are the approximate solutions to the CCFM-K equations. As a final remark, the parameter c = 1 in (91) leads to the replacement 1/b → c/b in all the formulae above.