Drell–Yan production with the CCFM-K evolution

Abstract

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\sim 1$ and transverse momentum ordering for $x \ll 1$. 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.

1 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 fixed-order 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\ll 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\sim M$, where fixed order perturbative QCD should apply. As it was shown in [15], when one considers moderate values of $M\sim 5-19\,\mathrm{GeV}$, the fixed-order predictions underestimate data in the region of $q_T\sim M$.

Such apparent troubles with the collinear factorization prompts us to approach the DY process using more general concepts like $k_T$-factorization [16,17,18,19,20,21] which allow to address the issue of intrinsic transverse momentum of partons. It should be mentioned, however, that the $k_T$-factorization approach is less popular since the higher QCD corrections are much harder to obtain than in the collinear framework. Also, unlike collinear factorization, it lacks formal proof. Despite these facts, many theoretical and phenomenological analyses of the DY process were made [22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37]. The theoretical improvement was accompanied by many experimental results from fixed–target [38,39,40,41] and collider [42,43,44,45,46] experiments.

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\sim 1$ to the region of small $x\ll 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\sim 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 economical 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\sim 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 off-shell 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.

2 CCFM approach

2.1 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

$$\begin{aligned} P_1={\textstyle {\frac{1}{2}}}\sqrt{S}(1,0,0,1)\,\quad \quad P_2={\textstyle {\frac{1}{2}}}\sqrt{S}(1,0,0,-1) \end{aligned}$$

(1)

in which the momenta in the branching process shown in Fig. 1 are given by

$$\begin{aligned} k_{i-1}= & {} x_{i-1}P_1+\bar{x}_{i-1}P_2+k_{(i-1) T} \,,\quad \quad \nonumber \\ k_i= & {} x_i P_1+\bar{x}_i P_2+k_{iT}. \end{aligned}$$

(2)

Notice that the momentum fraction are proportional to the plus/minus components, e.g.

$$\begin{aligned} x_i=\frac{k_i^+}{P_1^+} = \frac{k_{i}^0 +k_i^3}{\sqrt{S}} \,,\quad \quad \bar{x}_i=\frac{k_i^-}{P_2^-} = \frac{k_{i}^0 -k_i^3}{\sqrt{S}}. \end{aligned}$$

(3)

The emitted parton momentum can be found from the momentum conservation

$$\begin{aligned} p_i= & {} k_{(i-1) T}-k_{iT}= (x_{i-1}-x_i)P_1+(\bar{x}_{i-1}-\bar{x}_{i})P_2\nonumber \\&+(k_{(i-1) T}-k_{iT}), \end{aligned}$$

(4)

where $x_{i-1}>x_i$. Denoting the transverse component by $p_{iT}=k_{(i-1) T}-k_{iT} =(0,{\vec {p}}_{iT},0)$ and assuming that $p_i^2=0$, one can compute the minus component to find

$$\begin{aligned} p_i=(x_{i-1}-x_i)\,P_1 +\frac{{\vec {p}}_{iT}^{\,2}}{(x_{i-1}-x_{i})S}\,P_2+p_{iT} \end{aligned}$$

(5)

The rapidity of the emitted parton is given by

$$\begin{aligned} y_i=\frac{1}{2}\ln \left( \frac{p_i^+}{p_i^-}\right) = -\ln \frac{|{\vec {p}}_{iT}|}{x_{i-1}(1-z_i)\sqrt{S}} , \end{aligned}$$

(6)

where $z_i={x_i}/{x_{i-1}}<1$. In the massless case $y_i=-\ln \tan ({\theta _i}/{2})$ where $\theta _i$ is the emission angle with respect to the z axis defined by the collinear momenta $P_1$ and $P_2$, therefore

$$\begin{aligned} \tan \frac{\theta _i}{2} = \frac{|{\vec {p}}_{iT}|}{x_{i-1}(1-z_i)\sqrt{S}}\,. \end{aligned}$$

(7)

Fig. 1

Parton branching momenta

The CCFM branching scheme [50] is defined with the help of the rescaled transverse momentum of the emitted parton,

$$\begin{aligned} {\vec {q}}_{iT}=\frac{{\vec {p}}_{iT}}{1-z_i}. \end{aligned}$$

(8)

Thus, the transverse momentum conservation at the vertex i reads

$$\begin{aligned} {\vec {k}}_{(i-1)T} = {\vec {k}}_{iT} +(1-z_i)\, {\vec {q}}_{iT}. \end{aligned}$$

(9)

From (7) we obtain for the modulus

$$\begin{aligned} q_{iT}\equiv |{\vec {q}}_{iT}|=x_{i-1}\sqrt{S}\tan ({\theta _i}/{2}). \end{aligned}$$

(10)

The single loop approximation is defined by the condition

$$\begin{aligned} q_{iT} > q_{(i-1)T}\,, \end{aligned}$$

(11)

thus for $z_i\rightarrow 0$ we find the transverse momentum ordering of emitted partons

$$\begin{aligned} |{\vec {p}}_{iT}|> |{\vec {p}}_{(i-1)T}|. \end{aligned}$$

(12)

For finite $z\sim 1$, however, from (10) applied to (11) we have

$$\begin{aligned} x_{i-1}\tan \frac{\theta _i}{2} > x_{i-2}\tan \frac{\theta _{i-1}}{2} \end{aligned}$$

(13)

and for $z_{i-1}=(x_{i-1}/x_{i-2})\rightarrow 1$ we obtain angular ordering of parton emissions

$$\begin{aligned} \theta _i > \theta _{i-1}. \end{aligned}$$

(14)

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\rightarrow 0$ and angular ordering for $z\rightarrow 1$. The all loop approximation is defined by the condition

$$\begin{aligned} q_{iT} > z_{i-1}\,q_{(i-1)T} \end{aligned}$$

(15)

which is equivalent to the condition

$$\begin{aligned} \tan \frac{\theta _i}{2} > \tan \frac{\theta _{i-1}}{2}\,, \end{aligned}$$

(16)

giving the angular ordering condition (14). Thus, in the all loop approximation partons are emitted with angular ordering for any value of z.

2.2 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

$$\begin{aligned} f_g(x,{k_T},Q)&= f_g^0(x,{k_T},Q_0)\nonumber \\&\quad +\int \frac{d^2{\vec {q}}}{\pi q^2}\int _x^1\frac{dz}{z}\,\theta (Q-zq)\theta (q-Q_0)\,\frac{\alpha _s(q)}{2\pi } \Delta _S(Q,zq)\ \nonumber \\&\quad \times (2N_c)\left[ \frac{\Delta _{NS}({k_T},q,z)}{z}+\frac{\theta (1-z-Q_0/q)}{1-z}\right] \nonumber \\&\quad \times f_g\left( \frac{x}{z},|{\vec {k}}_T+(1-z){\vec {q}}|, q\right) \end{aligned}$$

(17)

which relates the gluon distribution at vertex i with the gluon distribution at vertex $(i-1)$. In the above, ${k_T}=|{\vec {k}}_T|$ and $q=|{\vec {q}}|$ are transverse momenta depicted in Fig. 1 and $\Delta _S$ is the Sudakov form factor given by

$$\begin{aligned}&\Delta _S(Q,zq)\nonumber \\&\quad =\exp \left\{ -\int _{(zq)^2}^{Q^2}\frac{dp^2}{p^2}\frac{\alpha _s(p^2)}{2\pi }\int _0^{1-Q^0/p}dz^\prime z^\prime P_{gg}(z^\prime )\right\} ,\nonumber \\ \end{aligned}$$

(18)

where $P_{gg}$ is the gluon-gluon splitting function (27), which resums virtual corrections for $z\rightarrow 1$. For $z\rightarrow 0$, the virtual emissions are resumed by the non-Sudakov form factor

$$\begin{aligned} \Delta _{NS}({k_T},q,z) =\exp \bigg \{ -\int _{z}^{1}dz^\prime \,\frac{2N_c}{z^\prime } \int ^{k^2_T}_{(z^\prime q)^2}\frac{dp^2}{p^2} \frac{\alpha _s(p^2)}{2\pi }\bigg \}\,.\nonumber \\ \end{aligned}$$

(19)

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,

$$\begin{aligned} zq=z\left( \frac{x}{z}\right) \sqrt{S}\tan \frac{\theta }{2}<Q\,. \end{aligned}$$

(20)

Therefore, for given x and $\sqrt{S}$, the hard scale determines the maximal emission angle

$$\begin{aligned} \theta _{\mathrm{max}} = 2\arctan \frac{Q}{x\sqrt{S}}. \end{aligned}$$

(21)

The second theta function in (17) imposes the condition

$$\begin{aligned} q>Q_0\gg \Lambda _{QCD} \end{aligned}$$

(22)

which assures that $\alpha _s(q)\ll 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

$$\begin{aligned} 0<z<(1- {Q_0}/{q)}\,, \end{aligned}$$

(23)

which allows to avoid singularity of $P_{gg}(z)$ at $z=1$. Non-perturbative effects are encoded in the initial condition, $f^0_g(x,{k_T},Q_0)$, imposed at a scale $Q_0$.

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].

2.3 CCFM-K evolution equations

The mixing between the transverse and longitudinal variables in the theta function $\theta (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]

$$\begin{aligned} \Theta (Q-zq)~~&\rightarrow ~~ \Theta (Q-q), \nonumber \\ \Delta _{NS}(k_T,q,z)~~&\rightarrow ~~ 1. \end{aligned}$$

(24)

Thus, for $z\rightarrow 0$, the angular ordering is replaced by the transverse momentum ordering while for $z\rightarrow 1$ the angular ordering is still valid. In addition, quark splittings, $q\rightarrow qg$, $\bar{q}\rightarrow \bar{q}g$ and $g\rightarrow q{\overline{q}}$, are taken into account which allow to introduce quark distributions $f_{i=1,\ldots ,2N_f}$ in addition to the gluon distribution $f_g$. In this way we obtain [63]

$$\begin{aligned}&f_{i}(x,{k_T},Q) \nonumber \\&\quad = f_i^0(x,{k_T}) \nonumber \\&\qquad + \int _0^1 \frac{dz}{z}\int \frac{d^2{\vec {q}}}{\pi q^2}\frac{\alpha _s(q^2)}{2\pi }\theta (Q-q)\theta (q-Q_0) \nonumber \\&\qquad \times \left\{ \theta (z-x)\left[ P_{qq}(z) f_i\left( \frac{x}{z},k_T^\prime ,q\right) \right. \right. \nonumber \\&\qquad \left. +P_{qg}(z)f_g\left( \frac{x}{z},k_T^\prime ,q\right) \right] \nonumber \\&\qquad \left. -z P_{qq}(z) f_i(x,{k_T},q) \right\} , \nonumber \\&f_g(x,{k_T},Q) \nonumber \\&\quad = f_g^0(x,{k_T}) + \int _0^1 \frac{dz}{z}\int \frac{d^2{\vec {q}}}{\pi q^2}\frac{\alpha _s(q^2)}{2\pi }\theta (Q-q)\theta (q-Q_0) \nonumber \\&\qquad \times \left\{ \theta (z-x)\left[ P_{gq}(z)\sum _{i=1}^{2N_f} f_i\left( \frac{x}{z},k_T^\prime ,q\right) \right. \right. \nonumber \\&\qquad \left. +P_{gg}(z)f_g\left( \frac{x}{z},k_T^\prime ,q\right) \right] \nonumber \\&\qquad \left. - z\left[ P_{gg}(z) +2N_fP_{qg}(z)\right] f_g(x,{k_T},q) \right\} , \end{aligned}$$

(25)

where the argument of the parton distributions on the r.h.s. equals

$$\begin{aligned} k_T^\prime =|{\vec {k}}_T+(1-z){\vec {q}}|. \end{aligned}$$

(26)

The one loop real emission splitting functions are given by¹

$$\begin{aligned} \nonumber P_{qq}(z)&= C_F\frac{1+z^2}{1-z}, \\ \nonumber P_{qg}(z)&=T_R\left\{ z^2+(1-z)^2\right\} , \\ \nonumber P_{gq}(z)&=C_F\frac{1+(1-z)^2}{z}, \\ P_{gg}(z)&=2C_A\left\{ \frac{z}{1-z}+\frac{1-z}{z} +z(1-z)\right\} , \end{aligned}$$

(27)

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 $\vec {k}_T$, the ordinary DGLAP equations are found for the collinear quark and gluon distributions,

$$\begin{aligned} q_i(x,Q)= & {} \int {d^2{\vec {k}}_T}\,f_{i}(x,{k_T},Q)\,,\quad \quad \nonumber \\ g(x,Q)= & {} \int {d^2{\vec {k}}_T}\,f_{g}(x,{k_T},Q)\,, \end{aligned}$$

(28)

Equation (25) can be written with the help of the Fourier transformation

$$\begin{aligned} \tilde{f}_{i,g}(x,{\vec {b}},Q) = \int {d^2{\vec {k}}_T}\,{e}^{i{\vec {k}}_T\cdot {\vec {b}}}f_{i,g}(x,{k_T},Q) \end{aligned}$$

(29)

which for ${\vec {b}}=0$ gives the PDFs (28). Since the parton distributions depend on ${k_T}=|\vec {k}_T|$, we perform the azimuthal angle integration with help of the relation

$$\begin{aligned} {e}^{i\vec {k}_T\cdot {\vec {b}}}=J_0({k_T}b)+2\sum _{n=1}^\infty i^nJ_n({k_T}b) \cos \phi \end{aligned}$$

(30)

and obtain the parton distributions which depend on $b=|{\vec {b}}|$,

$$\begin{aligned} \tilde{f}_{i,g}(x,b,Q) = \pi \int _0^\infty dk_T^2\,J_0({k_T}b)f_{i,g}(x,{k_T},Q)\,. \end{aligned}$$

(31)

Thus, taking the Fourier transform of both sides of Eq. (25), we find the evolution equations which are diagonal in b:

$$\begin{aligned} \frac{\partial \tilde{f}_i(x,b,Q)}{\partial \ln Q^2}&= \frac{\alpha _s(Q^2)}{2\pi }\int _0^1\frac{dz}{z} \Bigg \{ \theta (z-x) J_0((1-z)Qb) \nonumber \\&\quad \times \left[ P_{qq}(z) \tilde{f}_i\left( \frac{x}{z},b,Q\right) \right. \nonumber \\&\quad \left. + P_{qg}(z)\tilde{f}_g\left( \frac{x}{z},b,Q\right) \right] \nonumber \\&\quad - zP_{qq}(z) \tilde{f}_i(x,b,Q) \Bigg \}, \nonumber \\ \frac{\partial \tilde{f}_g(x,b,Q)}{\partial \ln Q^2}&= \frac{\alpha _s(Q^2)}{2\pi }\int _0^1\frac{dz}{z} \Bigg \{ \theta (z-x) J_0((1-z)Qb) \nonumber \\&\quad \times \left[ P_{gq}(z) \sum _{i=1}^{2N_f}\tilde{f}_i\left( \frac{x}{z},b,Q\right) \right. \nonumber \\&\quad \left. + P_{gg}(z)\tilde{f}_g\left( \frac{x}{z},b,Q\right) \right] \nonumber \\&\quad - z\left[ P_{gg}(z) +2N_fP_{qg}(z)\right] \tilde{f}_g(x,b,Q) \Bigg \}. \end{aligned}$$

(32)

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.

$$\begin{aligned} \tilde{f}_{i}(x,b=0,Q)= q_i(x,Q),\quad \tilde{f}_{g}(x,b=0,Q)= g(x,Q)\nonumber \\ \end{aligned}$$

(33)

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 $\alpha _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 $\alpha _s$ is inside the integrals over z since it depends on the transverse momentum of an emitted parton, ${k_T}=(1-z) Q$. In such a case, a cutoff on the upper limit of z is necessary to avoid the Landau pole in $\alpha _s({k_T})$.

2.4 Initial conditions and b-dependence

Fig. 2

The singlet quark distribution $f_S^{CCFMK}$ from Eq. (36) as a function of b for fixed $x=10^{-2}$ and two scales: initial $Q_0=1~\mathrm{GeV}$ and final $Q=100~\mathrm{GeV}$ (left plot). These distributions where multiplied by the gaussian form factor (35) on the right plot

In order to solve Eq. (32), we need initial conditions specified as functions of x and b at some perturbative scale $Q_0\gg \Lambda _{QCD}$. 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

$$\begin{aligned} \tilde{f}_{i}(x,b,Q_0)= & {} q_i(x,Q_0)F(b)\,,\quad \quad \nonumber \\ \tilde{f}_{g}(x,b,Q_0)= & {} g(x,Q_0)F(b), \end{aligned}$$

(34)

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$,

$$\begin{aligned} F(b) = \exp (-b^2/b_0^2). \end{aligned}$$

(35)

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

$$\begin{aligned} f_{i,g}(x,b,Q)= f_{i,g}^{\mathrm{CCFMK}}(x,b,Q)\,F(b)\,. \end{aligned}$$

(36)

where $f_{i,g}^{\mathrm{CCFMK}}$ is the solution for $F(b)\equiv 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=\sum _{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~\mathrm{GeV}$ with the MSTW08 LO PDFs [68] while the solid curves are evolved to $Q=100~\mathrm{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].

2.5 Relation to CSS resummation

Fig. 3

The solution of Eq. (32) for singlet quark (left plot) and gluon distributions (right plot) at $Q=100~\mathrm{GeV}$ and $x=10^{-2}$, obtained for $F(b)\equiv 1$ (solid lines) and the CSS approximation (38) (dot-dashed lines). The ratio CSS/CCFMK curves are shown as the red solid lines

The CCFM-K equations contain a part of the Collins, Soper and Sterman (CSS) resummation of soft parton emissions in the limit $z\rightarrow 1$. In “Appendix A”, we present the proof that for the values of the parameter b such that

$$\begin{aligned} 1/Q\ll b\ll 1/Q_0\,, \end{aligned}$$

(37)

the solution to the CCFM-K equations (32) is given by the CSS formulas [7]

$$\begin{aligned} f_i(x,b,Q)&=\exp \bigg \{-\int _{1/b^2}^{Q^2}\frac{dq^2}{q^2}\frac{\alpha _s(q^2)}{2\pi }\nonumber \\&\quad \times \left[ A_{q}^{(1)}\ln (q^2b^2)+B_{q}^{(1)}\right] \bigg \}\,q_i(x,1/b) \nonumber \\ f_g(x,b,Q)&=\exp \bigg \{-\int _{1/b^2}^{Q^2}\frac{dq^2}{q^2}\frac{\alpha _s(q^2)}{2\pi }\nonumber \\&\quad \times \left[ A_{g}^{(1)}\ln (q^2b^2)+B_{g}^{(1)}\right] \bigg \}\,g(x,1/b) \end{aligned}$$

(38)

where the parameters $A_{q,g}^{(1)}$ and $B_{q,g}^{(1)}$ are defined in Eq. (90). The above formulas were derived by picking large logarithms, $\ln (Qb)$ and $\ln (1/Q_0b)$, in the CCFM solution. It should be noted, however, that Eq. (38) do not contain the NLL (next-to-leading logarithmic) terms proportional to $\alpha _s^2$ 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~\mathrm{GeV}$ and $Q=100~\mathrm{GeV}$, the range (37) corresponds to $b\in [10^{-2},1]$ GeV$^{-1}$. We see that the CSS approximations extracted from the CCFM-K equation works reasonable well for $b\in [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~\mathrm{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_{\mathrm{min}}\sim 1/Q$. Analyzing this idea in the context of the CCFM-K approach, however, would go beyond the main thrust of our analysis.

3 Drell–Yan cross section with ${k_T}$-dependent PDFs

The Drell–Yan cross section differential in photon’s momentum is given by

$$\begin{aligned} \frac{d\sigma ^{DY}}{dy_\gamma \,dM^2 \,d^2q_{T}} =\frac{\alpha _{\mathrm{em}}^2}{24\pi ^3S^2M^2} \left( -W^{\mu }_{~\mu }\right) , \end{aligned}$$

(39)

where $(y_\gamma , M,q_{T})$ are photon’s rapidity, virtuality and transverse momentum while $W^{\mu }_{~\mu }$ is the trace of the hadronic tensor $W^{\mu \nu }$. With the lowest order matrix element for the process $q{\overline{q}}\rightarrow \gamma ^*$, the trace is given by the transverse momentum factorization formula

(40)

where $\vec {k}_{1T},\vec {k}_{2T}$ are quark transverse momenta, $x_{1,2}$ are their longitudinal momenta and $f_{i}, {f}_{\bar{i}}$ are transverse momentum dependent quark/antiquark distributions taken at the scale $Q=M$.

3.1 On-shell matrix element

The trace in (40) is the squared matrix element of the process $q(k_1){\overline{q}}(k_2)\rightarrow \gamma ^*$ in the lowest order. For the on-shell matrix element, we use the quark/antiquark momenta in the collinear approximation

$$\begin{aligned} k_1=x_1P_1\,,\quad k_2=x_2P_2\,. \end{aligned}$$

(41)

what assures gauge invariance of the matrix elements. In such a case

(42)

and the DY cross section (39) is given by

$$\begin{aligned} \frac{d\sigma ^{DY}}{dy_\gamma \,dM^2 \,d^2q_{T}}&= \frac{4\pi \alpha _{\mathrm{em}}^2}{3N_c\,M^4} \int d^2k_{1T}d^2k_{2T}\,\delta ^2(\vec {q}_T-\vec {k}_{1T}-\vec {k}_{2T}) \nonumber \\&\quad \times \sum _{i=1}^{N_f} e_i^2\,x_1x_2 \left[ f_{i}(x_1,\vec {k}_{1T}, M)\, f_{\bar{i}}(x_2,\vec {k}_{2T},M)+(1\leftrightarrow 2)\right] , \end{aligned}$$

(43)

It is easy to check that after integrating (43) over $\vec {q}_T$, we find the leading order form of the Drell–Yan cross section with collinear PDFs given by Eq. (28)

$$\begin{aligned}&\frac{d\sigma ^{DY}}{dy_\gamma \,dM^2}\nonumber \\&\quad = \frac{\sigma _0}{M^4}\sum _{i=1}^{N_f} e_i^2\,x_1x_2\left[ q_i(x_1,M^2) \bar{q}_{{i}}(x_2,M^2)+(1\leftrightarrow 2)\right] ,\nonumber \\ \end{aligned}$$

(44)

where $\sigma _0 = {4\pi \alpha _{\mathrm{em}}^2}/{3N_c }$. Inserting the delta function

$$\begin{aligned} \delta ^2(\vec {k}_{1T}+\vec {k}_{2T}-\vec {q}_T)=\int \frac{d^2b}{(2\pi )^2}\,{e}^{i(\vec {k}_{1T}+\vec {k}_{2T}-\vec {q}_T)\cdot {\vec {b}}} \end{aligned}$$

(45)

to Eq. (43), we find the DY cross section with the b-dependent parton distributions (29)

$$\begin{aligned} \frac{d\sigma ^{DY}_{\mathrm{on-shell}} }{dy_\gamma \,dM^2 \,d^2q_{T}}&= \frac{\sigma _0}{M^4} \int \frac{d^2{\vec {b}}}{(2\pi )^2}\,{e}^{-i\vec {q}_T\cdot {\vec {b}}}\, \sum _{i=1}^{N_f} e_i^2\,x_1x_2\nonumber \\&\quad \times \left[ \tilde{f}_{i}(x_1,{\vec {b}}, M)\, {\tilde{f}}_{\bar{i}}(x_2, {\vec {b}},M)+(1\leftrightarrow 2)\right] . \end{aligned}$$

(46)

For the parton distributions which depend on $b=|{\vec {b}}|$, the angular integration with the help of relation (30) gives

$$\begin{aligned} \frac{d\sigma ^{DY}_{\mathrm{on-shell}} }{dy_\gamma \,dM^2 \,dq_{T}^2}&= \frac{\sigma _0}{M^4} \int _0^\infty \frac{bdb}{2}J_0(q_{T}b)\sum _{i=1}^{N_f} e_i^2\,x_1x_2\nonumber \\&\quad \times \left[ \tilde{f}_{i}(x_1, b, M)\, {\tilde{f}}_{{\bar{i}}}(x_2, b,M)+(1\leftrightarrow 2)\right] . \end{aligned}$$

(47)

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 on-shell form

$$\begin{aligned} x_{1,2}={\frac{M}{\sqrt{S}}}\,{e}^{\pm y_\gamma }\,. \end{aligned}$$

(48)

3.2 Off-shell matrix elements

In approach with the off-shell matrix, the trace (42) is replaced by

(49)

where $\Gamma ^\mu$ is the Fadin-Sherman photon-quark vertex [70, 71]

(50)

and the quark/antiquark momenta $k_{1,2}$ take into account transverse components

$$\begin{aligned} k_1=x_1P_1+k_{1T}\,,\quad k_2=x_2P_2+k_{2T}, \end{aligned}$$

(51)

They are given by ${k_T}_{i}=(0,{\vec {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

$$\begin{aligned} x_{1,2}={\frac{M_T}{\sqrt{S}}}\,{e}^{\pm y_\gamma }\,,\quad M_T=\sqrt{M^2+q_{T}^2}. \end{aligned}$$

(52)

It is easy to check that the Fadin-Sherman vertex obeys the gauge invariance relation

$$\begin{aligned} (k_1+k_2)_\mu \Gamma ^\mu (k_1,k_2)=0. \end{aligned}$$

(53)

Computing the trace (49), we obtain

(54)

where we used momentum conservation at the photon vertex to write the last equality. Notice that because of the transverse mass $M_\perp$ in the denominator, the off-shell kinematics takes into account the corrections in powers of $q^2_\perp /M^2$ to all orders. With thes results, the cross section (39) is given by

$$\begin{aligned} \frac{d\sigma ^{DY}_{\mathrm{off-shell}}}{dy_\gamma \,dM^2 \,d^2q_{T}}&= \frac{4\pi \alpha _{\mathrm{em}}^2}{3N_c\,M^2M_T^2} \nonumber \\&\quad \times \int d^2k_{1T}d^2k_{2T}\,\delta ^2(\vec {q}_T-\vec {k}_{1T}-\vec {k}_{2T}) \left( 1+\frac{\vec {k}_{1T}^2+\vec {k}_{2T}^2}{M^2}\right) \nonumber \\&\quad \times \sum _{i=1}^{N_f} e_i^2 x_1x_2\left[ f_{i}(x_1,\vec {k}_{1T}, M)\, f_{{\bar{i}}}(x_2,\vec {k}_{2T},M)+(1\leftrightarrow 2)\right] . \end{aligned}$$

(55)

Inserting the delta function (45) and performing the Fourier transformation, we obtain the following cross section with the parton distributions which depend on $b=|{\vec {b}}|$

$$\begin{aligned} \frac{d\sigma ^{DY}_{\mathrm{off-shell}}}{dy_\gamma \,dM^2 \,dq_{T}^2}&= \frac{\sigma _0}{M^2M_\perp ^2} \int _0^\infty \frac{bdb}{2}\,J_0(q_{T}b)\,\nonumber \\&\quad \times \sum _{i=1}^{N_f} e_i^2\,x_1x_2 \, \bigg \{\tilde{f}_{i}(x_1,b,M) \tilde{f}_{{\bar{i}}}(x_2,b,M) \nonumber \\&\quad -\frac{1}{M^2}\left( \Delta _b\tilde{f}_{i}(x_1,b,M)\, \tilde{f}_{{\bar{i}}}(x_2,b,M) \right. \nonumber \\&\quad \left. + f_{i}(x_1,b,M)\, \Delta _b \tilde{f}_{{\bar{i}}}(x_2,b,M)\right) +(1\leftrightarrow 2) \bigg \} \end{aligned}$$

(56)

where $\Delta _b$ is the radial part of the two-dimensional Laplacian

$$\begin{aligned} \Delta _b=\frac{\partial ^2}{\partial b^2}+\frac{1}{b} \frac{\partial }{\partial b}\,. \end{aligned}$$

(57)

By the comparison with the cross section (47), we see that (56) has different mass dependence,

$$\begin{aligned} d\sigma ^{DY}_{\mathrm{off-shell}} \sim \frac{\sigma _0}{M^4(1+q_{T}^2/M^2)},\quad \mathrm{vs.}\quad d\sigma ^{DY}_{\mathrm{on-shell}} \sim \frac{\sigma _0}{M^4}\ . \end{aligned}$$

(58)

It should be emphasized that the corrections $q_{T}^2/M^2$ which are resummed to the factor $1/M_\perp ^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.

Fig. 4

Transverse momentum dependence of the DY cross sections: data from proton–proton R209 experiment are compared with CCFM-K on-shell cross section (47) (solid line), CCFM-K off-shell cross section (56) (dash-dotted line) and CSS cross section (59) (dashed line)

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\sigma ^{DY}_{\mathrm{off-shell}}$ is suppressed by a factor $M^2/M_\perp ^2$ in comparison to $d\sigma ^{DY}_{\mathrm{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\sigma ^{DY}_{\mathrm{off-shell}}$ at large $q_\perp$. 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~\mathrm{GeV}$ with the MSTW08 LO PDFs [68] and the form factor (35) with $b_0=2.7~\mathrm{GeV}^{-1}$. The dash-dotted 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.

3.3 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}\rightarrow 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

$$\begin{aligned} {d \sigma ^{DY}_{CSS} \over d y_\gamma dM^2 d q_T^2}= & {} \frac{\sigma _0}{M^2} \int _0^\infty \frac{bdb}{2}J_0(q_T b)\nonumber \\&\times \sum _{i=1}^{N_f} e_i^2\left\{ W_{i{\bar{i}}}(x_1,x_2,b,Q)+ (1\leftrightarrow 2)\right\} ,~~~~ \nonumber \\ \end{aligned}$$

(59)

where $Q=M$ and $x_{1,2}={M}{e}^{\pm y_\gamma }/{\sqrt{S}}$. The parton luminosities $W_{i{\bar{i}}}$ in the above read

$$\begin{aligned} W_{i{\bar{i}}}(x_1,x_2,b,Q)= & {} f_{i}^\prime (x_1,c/b_{*})\,f_{{\bar{i}}}^\prime (x_2,c/b_{*})\,e^{2S(b_*,Q)}\nonumber \\&\times \, W_{NP}(x_1, x_2,b,Q), \end{aligned}$$

(60)

where $f_{i/{\bar{i}}}^\prime$ are effective quark/antiquark distributions

$$\begin{aligned} f^\prime _{i/{\bar{i}}}(x,\mu )= & {} \int _x^1\frac{dz}{z} \left\{ C_q({x}/{z},\alpha _s(\mu )))\,q_{i/{\bar{i}}}(z,\mu )\right. \nonumber \\&\left. +C_g({x}/{z},\alpha _s(\mu ))\,g(z,\mu )\right\} \end{aligned}$$

(61)

with the $\overline{\mathrm{MS}}$ NLO collinear PDFs $q_{i/{\bar{i}}}(z)$ and g(z) and the coefficient functions

$$\begin{aligned} C_{q}(z,\alpha _s)&=\delta (1-z)\nonumber \\&\quad +\frac{\alpha _s}{2\pi }C_F\left[ 1-z+\left( \frac{\pi ^2}{2}-4\right) \delta (1-z)\right\} , \nonumber \\ C_{g}(z,\alpha _s)&= \frac{\alpha _s}{2\pi }T_R\left[ 2z(1-z)\right] . \end{aligned}$$

(62)

The scale $\mu$ in (61) is given by $\mu =c/b_*$ with $c=2{e}^{-\gamma _E}\approx 1.12$ and

$$\begin{aligned} b_*=\frac{b}{\sqrt{1+{b^2}/{b_{\mathrm{max}}^2}}}. \end{aligned}$$

(63)

TIn this way, $b_*$ interpolates between $b_*=0$ and $b_*=b_{\mathrm{max}}$ for $b\rightarrow \infty$ such that the scale

$$\begin{aligned} \mu \in [c/b_{\mathrm{max}},\infty ) \end{aligned}$$

(64)

Thus, choosing $b_{\mathrm{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~\mathrm{GeV}$.

The power S in the exponent in (60) is given by

$$\begin{aligned} S(b,Q)\,= & {} \,-\int _{c^2/b^2}^{Q^2}{d q^2 \over q^2}\left[ A_q(\alpha _s(q^2))\ln \left( {Q^2\over q^2}\right) \right. \nonumber \\&\left. + B_q(\alpha _s(q^2))\right] \end{aligned}$$

(65)

where the coefficients $A_q$ and $B_q$ are defined by the general perturbative expansion

$$\begin{aligned} A_q(\alpha _s)= & {} \sum _{n=1}^{\infty }\left( {\alpha _s\over 2\pi }\right) ^nA_q^{(n)} \,,~~~~~~~~~~~\nonumber \\ B_q(\alpha _s)= & {} \sum _{n=1}^{\infty }\left( {\alpha _s\over 2\pi }\right) ^nB_q^{(n)}\,. \end{aligned}$$

(66)

Introducing $B=\ln (Q^2b^2/c^2)$ and $L=L(Q)=\ln (Q^2/\Lambda ^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]

$$\begin{aligned} A_{q}^{(1)}=C_F,\quad A_{q}^{(2)}= C_FK,\quad B_{q}^{(1)}=-{\frac{3}{2}}C_F \end{aligned}$$

(67)

and $K=C_A(\textstyle {\frac{67}{18}}-\textstyle {\frac{1}{6}}\pi ^2)-\textstyle {\frac{10}{9}}T_RN_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_q^{(2)}$, 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

$$\begin{aligned} \alpha _s(q^2)=\frac{1}{\beta _0 L(q)} -\frac{\beta _1}{\beta _0^3}\frac{\ln L(q)}{L^2(q)} \end{aligned}$$

(68)

where $\beta _0=(11C_A-4T_RN_f)/12\pi$ and $\beta _1=(153-19N_f)/24\pi ^2$, and performing the integration in (65), one obtains the final NLL form of S [74], which we use in our presentation

$$\begin{aligned} S&= \frac{A_q^{(1)}}{\pi \beta _0}\left[ L\ln \left( 1-\frac{B}{L}\right) +B\right] \nonumber \\&\quad - \frac{A_q^{(2)}}{\pi ^2\beta _0^2}\left[ \ln \left( 1-\frac{B}{L}\right) +\frac{B}{L-B}\right] +\frac{B_q^{(1)}}{\pi \beta _0} \ln \left( 1-\frac{B}{L}\right) \nonumber \\&\quad + \frac{A_q^{(1)}\beta _1}{\pi \beta _0^3}\left[ \frac{B}{L-B} (1+\ln L)\nonumber \right. \\&\left. \quad +\left( \frac{L}{L-B}+\ln L\right) \ln \left( 1-\frac{B}{L}\right) +\frac{1}{2} \ln ^2\left( 1-\frac{B}{L}\right) \right] . \end{aligned}$$

(69)

The factor $W_{NP}=\exp \{-S_{NP}\}$ in (60) describes the non-perturbative contribution [75,76,77]. In our presentation, we use the form factor from the BLNY fit [78] to the DY data:

$$\begin{aligned} S_{NP} = [a_1+a_2\ln ({Q}/{Q_1}) + a_3\ln (100\, x_1x_2)]\, b^2, \end{aligned}$$

(70)

where

$$\begin{aligned} a_1= & {} 0.21~\mathrm{GeV}^{2}\,,\quad a_2=0.68~\mathrm{GeV}^{2},\nonumber \\ a_3= & {} -0.1~\mathrm{GeV}^{2}\,,\quad Q_1=3.2~\mathrm{GeV}. \end{aligned}$$

(71)

4 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~\mathrm{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_{\mathrm{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\,\mathrm{GeV}$ makes the curves stronger suppressed for large $q_T$.

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

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

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_T$ – factorization) 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\sim M$).

4.1 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 $E d^3 \sigma /d^3 q$ measured in these experiments is related to (47) and (59) as follows

$$\begin{aligned} E \frac{d^3 \sigma }{d^3 q} =\frac{2 M\Delta M}{\pi } \frac{d\sigma ^{DY}}{dy_\gamma \,dM^2 \,dq_{T}^2} \end{aligned}$$

(72)

where $\Delta 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

$$\begin{aligned} x_F \equiv \frac{2q^3}{\sqrt{S}}= \frac{2\sqrt{M^2+q^2_T}}{\sqrt{S}} \sinh {y_\gamma }, \end{aligned}$$

(73)

where the r.h.s. gives the relation between $x_F$ and the DY photon rapidity $y_\gamma$. 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 $\sqrt{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$)². 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 $\sqrt{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 $\sqrt{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 $\sqrt{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 off-shell (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 we fixed just one free parameter $b_0$ in the form factor (35) to $b_0=2.7\,\mathrm{GeV}^{-2}$.

Fig. 9

Transverse momentum dependence of DY cross-section in proton–proton collisions: data from PHENIX [46] compared with on-shell CCFM-K (blue solid), off-shell CCFM-K (blue dash-dotted) and CSS (red dashed) approaches

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 \sim 9\,\mathrm{GeV}$, the data from E288 are significantly above the theoretical predictions which is related to the production of $\Upsilon$ 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 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\sim M$, as one should expect, see discussion at the end of Sect. 3.2.

4.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 $\sqrt{S}= 62$ GeV and PHENIX [46] with $\sqrt{S}= 200$ GeV. For R209 we apply a change of variables,

$$\begin{aligned} \frac{d \sigma }{d^2 q_T}= \int \limits _{5-8\text { GeV}} d M \frac{M\sqrt{S}}{\sqrt{M^2 + q_T^2}} \frac{d\sigma ^{DY}}{dy\,dM^2 \,dq_T^2}, \end{aligned}$$

(74)

whereas PHENIX is using the cross section $E d^3 \sigma /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\sim 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.

5 Summary

Using the CCFM-K parton distributions and the partonic cross-section 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.

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

$$\begin{aligned} \bar{f}(n,b,Q)=\int _0^1 dx\,x^{n} \tilde{f}(x,b,Q) \end{aligned}$$

(75)

to both sides of Eq. (32) to obtain

$$\begin{aligned} \frac{\partial \bar{f}_i(n,b,Q)}{\partial \ln Q^2}&= \bar{P}_{qq}(n,b,Q)\,\bar{f}_i(n,b,Q)\nonumber \\&\quad + \bar{P}_{qg}(n,b,Q)\,\bar{f}_g (n,b,Q), \nonumber \\ \frac{\partial \bar{f}_g(n,b,Q)}{\partial \ln Q^2}&= \bar{P}_{gq}(n,b,Q)\, \bar{f}_S(n,b,Q) \nonumber \\&\quad + \bar{P}_{gg}(n,b,Q)\, \bar{f}_g(n,b,Q), \end{aligned}$$

(76)

where $\bar{f}_S=\sum _{i=1}^{2N_f} \bar{f}_i$ is the singlet quark distribution and the Mellin moments of the splitting functions read

$$\begin{aligned} \bar{P}_{qq}(n,b,Q)&= \frac{\alpha _s(Q^2)}{2\pi } \int _0^1dz\,P_{qq}(z)\nonumber \\&\quad \times \left\{ z^n J_0((1-z)Qb)-1\right\} , \nonumber \\ \bar{P}_{qg} (n,b,Q)&= \frac{\alpha _s(Q^2)}{2\pi } \int _0^1dz\,P_{qg}(z)\,z^n J_0((1-z)Qb), \nonumber \\ \bar{P}_{gq}(n,b,Q)&= \frac{\alpha _s(Q^2)}{2\pi }\int _0^1dz\,P_{gq}(z) \,z^n J_0((1-z)Qb), \nonumber \\ \bar{P}_{gg} (n,b,Q)&= \frac{\alpha _s(Q^2)}{2\pi } \int _0^1dz \left\{ P_{gg}(z)\, z^n J_0((1-z)Qb)\nonumber \right. \\&\quad \left. -z[P_{gg}(z)+2N_fP_{qg}(z)]\right\} . \end{aligned}$$

(77)

Notice that for $b=0$ in (76) we obtain the DGLAP evolution equations. This fact motivates the following decomposition of the diagonal splitting functions

$$\begin{aligned} \bar{P}_{qq}(n,b,Q)&= \bar{P}_{qq1}(n,b,Q) + \bar{P}_{qq2}(b,Q), \nonumber \\ \bar{P}_{gg}(n,b,Q)&= \bar{P}_{gg1}(n,b,Q) + \bar{P}_{gg2}(b,Q), \end{aligned}$$

(78)

where

$$\begin{aligned} \bar{P}_{qq1}(n,b,Q)&= \frac{\alpha _s(Q^2)}{2\pi } \int _0^1dz\,P_{qq}(z)(z^n-1) J_0((1-z)Qb), \end{aligned}$$

(79)

$$\begin{aligned} \bar{P}_{qq2}(b,Q)&= \frac{\alpha _s(Q^2)}{2\pi } \int _0^1dz\,P_{qq}(z)\left\{ J_0((1-z)Qb)-1\right\} , \end{aligned}$$

(80)

$$\begin{aligned} \bar{P}_{gg1}(n,b,Q)&= \frac{\alpha _s(Q^2)}{2\pi } \int _0^1dz \left\{ z^n P_{gg}(z)- z[P_{gg}(z) \nonumber \right. \\&\quad \left. +2N_fP_{qg}(z)]\right\} J_0((1-z)Qb), \end{aligned}$$

(81)

$$\begin{aligned} \bar{P}_{gg2}(b,Q)&= \frac{\alpha _s(Q^2)}{2\pi } \int _0^1dz \,z[P_{gg}(z) + 2N_f P_{qg}(z)]\nonumber \\&\quad \times \left\{ J_0((1-z)Qb) -1\right\} . \end{aligned}$$

(82)

Thus, for $b=0$, $\bar{P}_{qq2}=\bar{P}_{gg2}=0$, and $\bar{P}_{qq1}$ and $\bar{P}_{gg1}$ become the ordinary Altarelli-Parisi splitting functions. This is why we write (76) in the form

$$\begin{aligned} \nonumber \frac{\partial \bar{f}_i}{\partial \ln Q^2}- \bar{P}_{qq2}\, \bar{f}_i&= \bar{P}_{qq1}\,\bar{f}_i + \bar{P}_{qg}\,\bar{f}_g, \\ \frac{\partial \bar{f}_g}{\partial \ln Q^2}- \bar{P}_{gg2}\, \bar{f}_g&= \bar{P}_{gq}\, \bar{f}_S + \bar{P}_{gg1}\, \bar{f}_g, \end{aligned}$$

(83)

where for simplicity we suppressed the arguments. We look for the solutions in the form

$$\begin{aligned} {\bar{f}}_i(n,b,Q)= & {} {e}^{S_q(b,Q)}\bar{\bar{f}}_i(n,b,Q)\,,~~~~~~~~~\nonumber \\ {\bar{f}}_g(n,b,Q)= & {} {e}^{S_g(b,Q)}\bar{\bar{f}}_g(n,b,Q), \end{aligned}$$

(84)

where

$$\begin{aligned} S_q(b,Q)= & {} \int _{Q_0^2}^{Q^2}\frac{dq^2}{q^2} \bar{P}_{qq2}(b,q),\nonumber \\ S_g(b,Q)= & {} \int _{Q_0^2}^{Q^2}\frac{dq^2}{q^2} \bar{P}_{gg2}(b,q). \end{aligned}$$

(85)

Inserting (84) to (83), we find

$$\begin{aligned} \frac{\partial \bar{\bar{f}}_i}{\partial \ln Q^2}&= \bar{P}_{qq1}\,\bar{\bar{f}}_i + \left( {e}^{S_g-S_q}\bar{P}_{qg}\right) \bar{\bar{f}}_g, \nonumber \\ \frac{\partial \bar{\bar{f}}_g}{\partial \ln Q^2}&= \left( {e}^{S_q-S_g}\bar{P}_{gq}\right) \bar{\bar{f}}_S +\bar{P}_{gg1}\, \bar{\bar{f}}_g. \end{aligned}$$

(86)

We are interested in the approximate solution when b obeys relation (37), i.e. for

$$\begin{aligned} Q_0\ll 1/b \ll Q\,. \end{aligned}$$

(87)

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)

$$\begin{aligned} \bar{f}_i(n,b,Q)&=\exp \bigg \{-\int _{1/b^2}^{Q^2}\frac{dq^2}{q^2}\frac{\alpha _s(q^2)}{2\pi }\nonumber \\&\quad \times \left[ A_{q}^{(1)}\ln (q^2b^2)+B_{q}^{(1)}\right] \bigg \}\,\bar{q}_i(n,1/b), \end{aligned}$$

(88)

$$\begin{aligned} \bar{f}_g(n,b,Q)&=\exp \bigg \{-\int _{1/b^2}^{Q^2}\frac{dq^2}{q^2}\frac{\alpha _s(q^2)}{2\pi }\nonumber \\&\quad \times \left[ A_{g}^{(1)}\ln (q^2b^2)+B_{g}^{(1)}\right] \bigg \}\,\bar{g}(n,1/b) \end{aligned}$$

(89)

where the coefficients

$$\begin{aligned} A_{q}^{(1)}= & {} C_F, \quad B_{q}^{(1)}=-\frac{3}{2}C_F,\nonumber \\ A_{g}^{(1)}= & {} C_A,\quad B_{g}^{(1)}=\frac{2}{3}T_RN_f-\frac{11}{6} C_A. \end{aligned}$$

(90)

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

$$\begin{aligned} J_0(u)\approx \theta (c-u),\quad 1- J_0(u) \approx \theta (u-c), \end{aligned}$$

(91)

where $\theta$ is the Heaviside step function and $c\sim 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) with $\bar{P}_{qq2}$ given by (80) is given by

$$\begin{aligned}&S_q(b,Q) = -\int _{Q_0^2}^{Q^2}\frac{dq^2}{q^2} \frac{\alpha _s(q^2)}{2\pi }\nonumber \\&\quad \int _0^1 dz\,P_{qq}(z)\,\theta [(1-z)q b-1]. \end{aligned}$$

(92)

For $qb<1$, the argument of the theta function is negative and $S_q=0$. Thus, we have

$$\begin{aligned} S_q(b,Q)= & {} - \int _{Q_0^2}^{Q^2}\frac{dq^2}{q^2} \frac{\alpha _s(q^2)}{2\pi } \theta (q -1/b)\nonumber \\&\times \int _0^{1} dz\,P_{qq}(z)\,\theta (1-z-1/qb). \end{aligned}$$

(93)

From the first theta function $q>1/b \gg Q_0$ which sets the lower integration limit to 1/b. In this way, we avoid resummation of large logarithms $\ln (1/Q_0b)$ which are shifted to the functions $\bar{\bar{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

$$\begin{aligned} P_{qq}(z)=C_F\left( \frac{2}{1-z}-(1+z)\right) , \end{aligned}$$

(94)

we find the result which agrees with the exponent in (88)

$$\begin{aligned} S_q(b,Q)&= - \int _{1/b^2}^{Q^2}\frac{dq^2}{q^2} \frac{\alpha _s(q^2)}{2\pi } \nonumber \\&\qquad \times \int _0^{1-1/q b} dz\,C_F\left\{ \frac{2}{1-z}-(1+z)\right\} \nonumber \\&\quad \approx -\int ^{Q^2}_{1/b^2}\frac{dq^2}{q^2} \frac{\alpha _s(q^2)}{2\pi } \left\{ C_F \ln (q^2 b^2) -\frac{3}{2}C_F\right\} , \end{aligned}$$

(95)

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)

$$\begin{aligned} S_g(b,q)\approx & {} -\int ^{Q^2}_{1/b^2}\frac{dq^2}{q^2} \frac{\alpha _s(q^2)}{2\pi } \left\{ C_A \ln (q^2 b^2)\right. \nonumber \\&\left. +\frac{2}{3}T_RN_f -\frac{11}{6}C_A \right\} . \end{aligned}$$

(96)

The large logs $\ln (1/Q_0b)$ are resummed using the DGLAP evolution equations. To show this, we write (86) in the integral form

$$\begin{aligned} \bar{\bar{f}}_i(n,b,Q)&= \bar{\bar{f}}_i(n,b,Q_0) \nonumber \\&\quad + \int _{Q_0^2}^{Q^2} \frac{d\mu ^2}{\mu ^2}\bar{P}_{qq1}(n,b,\mu )\,\bar{\bar{f}}_i(n,b,\mu ) \nonumber \\&\quad + \int _{Q_0^2}^{Q^2} \frac{d\mu ^2}{\mu ^2}{e}^{S_g(b,\mu )-S_q(b,\mu )} \bar{P}_{qg}(n,b,\mu )\,\bar{\bar{f}}_g(n,b,\mu )\,. \end{aligned}$$

(97)

The first integral on the r.h.s. is given by

$$\begin{aligned} I_1= & {} \int _{Q_0^2}^{Q^2} \frac{d\mu ^2}{\mu ^2} \frac{\alpha _s(\mu ^2)}{2\pi }\nonumber \\&\times \left( \int _0^1dz\,P_{qq}(z)(z^n-1) J_0((1-z)\mu b)\right) \bar{\bar{f}}_i(n,b,\mu )\nonumber \\ \end{aligned}$$

(98)

where we used (79). Approximating

$$\begin{aligned} J_0((1-z)\mu b)\approx \theta (1-(1-z)\mu b)\,, \end{aligned}$$

(99)

we notice that in order to get the leading logs $\ln (1/Q_0b)$ we have to assume that $\mu b< 1$. In such a case, the integration over z in (98) is not constrained but the integration over $\mu$ is limited to $\mu <1/b\ll Q$. In this way

$$\begin{aligned} I_1\approx \int _{Q_0^2}^{1/b^2} \frac{d\mu ^2}{\mu ^2}\frac{\alpha _s(\mu ^2)}{2\pi } \left( \int _0^1dz\,(z^n-1)P_{qq}(z)\right) \bar{\bar{f}}_i(n,b,\mu )\,.\nonumber \\ \end{aligned}$$

(100)

Applying the same approximations to second integral in (97), we obtain

$$\begin{aligned} I_2\approx & {} \int _{Q_0^2}^{1/b^2} \frac{d\mu ^2}{\mu ^2}\frac{\alpha _s(\mu ^2)}{2\pi }\,{e}^{S_q(b,\mu )-S_g(b,\mu )}\nonumber \\&\times \left( \int _0^1dzz^nP_{qg}(z)\right) \bar{\bar{f}}_g(n,b,\mu ) . \end{aligned}$$

(101)

From (95), in the integration region, $\mu <1/b$, we have $S_q(b,\mu )=S_g(b,\mu )=0$. Therefore, we can set the exponent in $I_2$ equal one and (97) reads

$$\begin{aligned} \bar{\bar{f}}_i(n,b,1/b)&= \bar{\bar{f}}_i(n,b,Q_0) + \int _{Q_0^2}^{1/b^2} \frac{d\mu ^2}{\mu ^2}\frac{\alpha _s(\mu ^2)}{2\pi } \nonumber \\&\qquad \times \left( \int _0^1dz\,P_{qq}(z)(z^n-1)\right) \,\bar{\bar{f}}_i(n,b,\mu ) ~+\nonumber \\&\qquad + \int _{Q_0^2}^{1/b^2} \frac{d\mu ^2}{\mu ^2}\frac{\alpha _s(\mu ^2)}{2\pi }\, \left( \int _0^1dzz^nP_{qg}(z)\right) \nonumber \\&\qquad \times \bar{\bar{f}}_g(n,b,\mu )\,.~~~~~~~ \end{aligned}$$

(102)

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 $\mu >1/b$, the theta function (99) imposes the condition $z>1-1/(\mu b)$ leads to subleading logarithms which we neglect. $\square$

To make a connection with the collinear PDFs, we note that for the values of b which we consider, $b\ll 1/Q_0\sim 1~\mathrm{GeV}^{-1}$, to good approximation $b\approx 0$ in the functions $\bar{\bar{f}}_{i,g}$ in (102). Thus, choosing the initial conditions equal to the Mellin moments of the collinear PDFs,

$$\begin{aligned} \bar{\bar{f}}_i(n,b\approx 0,Q_0)= & {} {\overline{q}}_i(n,Q_0),\nonumber \\ \bar{\bar{f}}_g(n,b\approx 0,Q_0)= & {} \bar{g}(n,Q_0), \end{aligned}$$

(103)

we find the collinear PDFs at the scale 1/b

$$\begin{aligned} \bar{\bar{f}}_i(n,b\approx 0,1/b)= & {} {\overline{q}}_i(n,1/b),\nonumber \\ \bar{\bar{f}}_g(n,b\approx 0,1/b)= & {} \bar{g}(n,1/b). \end{aligned}$$

(104)

This concludes the proof that (88) and (89) are the approximate solutions to the CCFM-K equations. As a final remark, the parameter $c\ne 1$ in (91) leads to the replacement $1/b\rightarrow c/b$ in all the formulae above.