The analysis of Drell-Yan lepton pair production in the P-P(\={P}) colliders using different angular ordering constraints and $k_t$-factorization approach

In this work, the P-P(\={P}) Drell-Yan lepton pair production (DY) differential cross sections at hadrons colliders, such as LHC and TEVATRON, are studied in the ${k_{t}}$-factorization framework. In order to take into account the transverse momenta of incoming partons, we use the unintegrated parton distribution functions (UPDF) of Kimber et al (KMR) and Martin et al (MRW) in the leading order (LO) and next-to-leading-order (NLO) levels with the input MMHT2014 PDF libraries. Based on the different off-shell partonic matrix elements, we analyze the behaviors of DY differential cross sections with respect to the invariant mass, the transverse momentum and the rapidity as well as the specific angular correlation between the produced leptons. The numerical results are compared with the experimental data, in different energies, which are reported by various collaborations, such as CDF, CMS, ATLAS and LHCb. It is shown that the NLO-MRW and KMR schemes predict closer results to the data compared to the LO-MRW, since we do not have fragmentation.


I. INTRODUCTION
Traditionally, the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equations [1-4] approach is used to obtain the quark, anti-quark and gluon densities, i.e., the parton distribution functions (PDF), a(x, µ 2 ). These functions depend on the Bjorken variable x and the hard scale µ 2 , and can be easily used in the collinear QCD factorization formalisms. In these DGLAP evolution approaches, the transverse momentum (k t ) components of partons are integrated over and there is not any degree of freedom for the initial gluon radiations and the transverse momentum (k t ) of the partons in the PDF. But, the outcome of hadron-hadron colliders at high energies indicates that the explicit inclusion of the intrinsic transverse momentum of initial hadron constituents is important to get accurate results and predictions. Therefore, the important inputs are the transverse momentum dependent parton distribution functions or the so-called "unintegrated" PDF (UPDF).
Theoretically, various methods are proposed to generate these fundamental quantities, i.e., UPDF, and among them the BFKL [5-9] (which is valid for the small x and the scale k 2 t ) and CCFM [10][11][12][13][14] (which is applicable at both the small and large x, the scales k 2 programs, such as ResBos [40], DYNNLO [41] and POWHEG+PYTHIA [42] event generators. The ResBos method simulates the vector-boson production and its decay, using a resumed treatment of the soft-gluon emissions at the NLO-logarithm (NNLL) and the γ * and Z/γ * contributions are simulated at the NLO accuracies. The DYNNLO approach is a parton level Monte Carlo program that computes the cross sections for the vector boson production in the P-P(P) collisions up to the NNLO in the pQCD theory. The PYTHIA program generates the LO QCD interactions via its parton shower algorithms. In the recent investigation [33] the distribution of dilepton transverse momentum and angular variable φ * η were calculated perturbatively at √ s = 8 TeV using the ResBos Monte Carlo generator at the NNLO accuracy and compared to the ATLAS data. Although, the results at low values of p T and φ * η show a good agreement with data, but this is not the case at high values of p T and φ * η . In the present report, it is intended to calculate the DY differential cross sections based on the KMR and MRW k t -factorization approaches using the corresponding off-shell transition amplitudes. We consider the three sub-processes namely, (1) q * +q * → γ * /Z → l + + l − and (2) q * + g * → γ * /Z → l + + l − + q and (3) q * +q * → γ * /Z + g → l + + l − + g at the LO and NLO levels, respectively. The dependence of the DY differential cross sections on the dilepton transverse momentum, the invariant mass and the rapidity distributions as well as the angular correlation between the produced leptons are calculated in the above frameworks and compared with the experimental data developed by the CDF, CMS, ATLAS and LHCb collaborations in both the Tevatron and LHC energies. The Harland-Lang et al.
(MMHT2014) PDF libraries [43] in the LO for both the KMR [15] and MRW [16] formalisms and the NLO for NLO-MRW approaches are considered. These calculations are performed for the off-shell incoming partons.
We should also point out here that, recently, the consideration of various angular ordering, as well as the generation of different UPDF becomes the subject of several reports [44,45].
In the reference [44], it is pointed out that the KMR UPDF, which are generated using the differential and integral approaches, give different results in the region where k t > µ.
Therefore it is concluded that, the integral form of the KMR UPDF gives the correct result, while for the application of differential form, one should use the cutoff dependent PDF. On the other hand, in the reference [46], the above idea is rejected, and a new term is added to the Sudakov form factor via a Heaviside step function, to set the Sudakov form factor equal to 1 in the region where k t > µ. Further, it is claim that the above two forms of KMR approach give the same result, i.e., there is no need to introduce cutoff dependent PDF.
Finally, in the reference [46], by referring to this report, i.e., [47] ( in which the predictions of the KMR approach with the AOC overestimates the data in case of the heavy quark production), it is suggested that the above problem is due to the freedom of parton to have transverse momentum larger than µ, and concluded that it is much suitable to use the KMR UPDF with the strong ordering constraint (SOC), which harshly cuts the transverse momentum in k t > µ region. Because of the above statements, and all of the problems appear in k t > µ, we compute the DY differential cross sections with respect to M ll and p ll T using the SOC KMR UPDF to check the sensitivity of our results in this region (see the figures 1 (panel f) and 2 (panel d)). One should also note that, as it is discussed in the reference [48], the result of k t -factorization should not be as good as collinear factorization approaches, in covering the experimental data, on the other hand, it is more simplistic, considering computer time consuming.
It should be also noted that the Sudakov form factor of the KMR approach does not obey the multiplication law according to the reference [49], but despite of this fact, it is interesting to point out that the normalization condition (see the equation (1)) is approximately satisfied in the KMR formalism, which is a critical issue in constructing any new UPDF.
The k t -factorization calculations were also performed by considering one or two of the above three sub-processes with the MSTW2008 PDF [50][51][52]. In these works, although the authors declare that they use the KMR formalism, but they do not take into account the factor 1/k 2 t in the cross section nor in the normalization formulas: Beside these, they use different angular ordering conditions with respect to the KMR prescriptions (we refer to them as semi-KMR). However, their results are surprisingly close to the experimental data. In the above references [51,52], it is claimed that the second and the third of the above sub-processes can be omitted by effectively using only the reggeized (off-shell) quark approach in the first sub-process. A brief discussion about the result of these reports and the comparison with our predictions are presented in the section III. In the reference [51], although the off-shell initial quarks are used, it is shown that utilizing the reggeized model and the effective vertexes guaranteed the gauge invariance of the transition matrix elements (TME). However, in our previous work [53], we showed that using the offshell initial quarks in the k t -factorization dynamics and in the small x regions leads to the gauge invariance of the TME, too.
To check the validity of our calculated cross sections the KaTie parton-level event generator [54] . The UPDF of this form are rarely investigated in the phenomenological applications of the k t -factorization. We include the sub-process q * +q * → γ * /Z → l + + l − + g which usually is neglected, e.g. [50]. In the others works, including those that are cited in our report [50], incorrectly, the combinations of KMR and MRW formalisms are used and they forget about the importance of the normalization constraint (the equation (1)) on the UPDF. This point is discussed in details in the reference [53]. One should note that the KMR prescription is a semi-NLO. Another important item is the fragmentation effect, which does not present in the processes that are discussed in our report. Because of that, as it is explained in the paper, the KMR and NLO-MRW procedures demonstrate better agreement to the experimental data.
The outline of our paper is as following. In the section II, we briefly present the basic cross section formulas of k t -factorization (II A) approach and the derivation of input UPDF (II B).
In the section III we present numerical calculations (III A), results presentations (III B) and discussions (III B). The section IV expresses our conclusions.

II. THE THEORETICAL FRAMEWORK OF DY
A. The k t -factorization cross section formulas Our DY differential cross section calculations are based on the k t -factorization in the KMR and MRW UPDF [15,16] approaches. Therefore in this section, we describe the theoretical framework of these approaches as well as the corresponding matrix elements (also see the appendix A). As we pointed out in the introduction, we include all the subprocesses contributions up to the αα s levels, namely: q * +q * → γ * /Z → l + + l − , q * + g * → γ * /Z → l + + l − + q and q * +q * → γ * /Z + g → l + + l − + g. From the kinematical point of view, if we show the four-momenta of the incoming protons (partons) by P (1) (k (1) ) and P (2) (k (2) ) and neglect their masses, then in the proton center of mass framework we have: where √ s is the total center of mass energy. In the high energy and the leading-logapproximation kinematics, the corresponding partons four-momenta can be written in terms of their transverse momenta k 1t and k 2t and the fraction (x i ) of the incoming protons momentum as: There are some relations for the above three sub-processes due to the energy-momentum conservation law as following: x 1 = 1 √ s (m 1t e y 1 + m 2t e y 2 ), (5) for the first sub-process and x 1 = 1 √ s (m 1t e y 1 + m 2t e y 2 + m 3t e y 3 ), (8) for the second and third sub-processes, where p it , y i and m it are the transverse momenta, the rapidities and the transverse masses (m 2 it = m 2 i + p 2 it ) of the produced particles, (i=1 and 2 for leptons and i=3 for (anti-)quark or gluon), respectively.
To calculate the matrix elements squared in the k t -factorization framework [55], the summation over the incoming off-shell gluon polarizations is carried out as: where k t is the gluon transverse momentum. For the off-shell quarks spinors with momentum k (after imposing the Sudakov decomposition in the high energy and the leading-logapproximation kinematics [56]), we have u(k)ū(k) xP , where x represents the fractional longitudinal momentum of proton, see the references [53,57,58]. Also, the effective vertices are used to calculate the Feynman amplitudes to test and ensure the gauge invariance of the different matrix elements [59][60][61]. It is worth to point out that a similar technique is also developed, using the Slavnov-Taylor identities by the means of the helicity amplitude [62][63][64] and being checked against those obtained by usage of Lipatov's effective action [59][60][61], as we pointed out.
To calculate the differential cross sections of DY, according to the k t -factorization theorem, for 2 → 2 sub-process we have: k 2 2t dp 2 1t dp 2 2t dy 1 dy 2 dφ 1 2π and for 2 → 3 sub-process one finds: k 2 2t dp 2 1t dp 2 2t dy 1 dy 2 dy 3 dφ 1 2π where f q (x i , k 2 it , µ 2 ) are the UPDF, which depend on the two hard scales, k 2 t and µ 2 , and they can be written in terms of the usual PDF. As we pointed out in the introduction, in the present calculations, the MMHT2014 PDF [43] is used for calculating the UPDF. In the above formula, M j i are the off-shell matrix elements which are presented for the three different sub-processes in the appendix A. Note that when we squared the matrix element of each three sub-processes we get the interference effect between γ * and Z production, which will be discussed in the section III. The azimuthal angles of the initial partons and the produced leptons are presented by φ 1 and φ 2 , and ψ 1 and ψ 2 , respectively. Then the total cross section can be written as: To calculate the UPDF of (anti-)quarks and gluons in a proton, we apply the LO KMR, LO MRW and NLO-MRW approaches [15,16]. In the following each of them will be described.
In the KMR method the UPDF of each parton, which means the probability to find a parton with transverse momentum k t and fractional momentum x at hard scale µ 2 are given by: where which is the familiar Sudakov survival form factor and limits the emissions of partons between k 2 t and µ 2 scales [15,16]. P (0) ab (z) are the usual LO splitting functions. In this formula the angular-ordering constraint (AOC) [10-14, 65, 66], ∆, is applied in the upper limit of the integration, which is an infrared cutoff to prevent the soft gluon singularities arise from the splitting functions and defined as: Note that this constraint is imposed on both quark and gluon radiations. b x z , k 2 t are the LO PDF, and in this work they are taken from the MMHT2014 libraries [43].
To determine the UPDF we also apply the MRW prescription which is similar to the KMR formalism, but the AOC only acts on the terms which include the on shell gluon emissions. For the quarks and the gluons they take the following forms: with and with respectively. In the above equations, z max = 1 − z min = µ µ+kt [67]. By expanding MRW to the NLO level, we have: In this formalism the Sudakov form factor is defined as: The higher order splitting functions are presented in the appendix B.

A. Numerical calculations
In this section, we present the kinematics and theoretical inputs of our calculations.
First, we calculate the UPDF based on the different k t -factorization schemes by using two different methods, i.e., KMR and MRW. Through our calculations, we set the renormalization and factorization scales to be equal to is the invariant mass of produced dilepton and as usual we consider the default value ζ = 1 [53]. We let this parameter to vary from 1/2 to 2, to estimate the scale uncertainties of our calculations. We also set m Z = 91.187 GeV and Λ QCD = 200M eV with n f = 4 active quark flavors. Using the LO coupling constant, we get α s (M 2 Z ) = 0.123 (g W = 0.66). Second, with the massless quarks approximation, the calculation of transition matrix elements squared is carried out, using the small x approximation presented in the appendix A, by the means of FeynCalc [68], i.e., the mathematica package for symbolic semi-automatic evaluation of Feynman diagrams. In the present report, the non-logarithmic loop corrections to the q-q annihilation cross section are taken into account by applying the effective K-factor with a particular scale choice of µ 2 = p 4/3 T M 2/3 as it was done, for example, in the references [26,50,69], i.e., where p T , (p T =| p 1t + p 2t |, see above the equation (9)), is the transverse momentum of produced dilepton and C F is the color factor. To calculate the multidimensional integration, the VEGAS routine [70] is used. The differential cross sections at several center of mass energies, i.e., 1.96, 7 and 8 TeV as a function of the dilepton invariant mass (M ), rapidity (y), transverse momentum (p T ) and the variable φ * η [34,[71][72][73], i.e., are calculated, with φ acop = π−|∆φ|, where ∆η and ∆φ are the pseudorapidity and azimuthal angles differences between the produced leptons, as well, respectively. The variable φ * η is correlated to the quantity |p T |/M and both of them probes the same physics as the dilepton transverse momentum, but it gives a better experimental resolution [74][75][76].

B. Results presentations
The results of above numerical calculations are compared with the experimental data of DY at the Tevatron and LHC laboratories with the total center of mass energy √ s = 1.8 TeV and √ s = 7 and 8 TeV, respectively. We use the data from different groups such as the CDF, CMS, ATLAS and LHCb collaborations. The available pQCD predictions are also presented in each figure.
The above comparisons are demonstrated in the figures 1 to 10 as follows: (1): In all of the figures, the numerical results related to the KMR UPDF are shown in the left panels in which the dash, dotted-dash and dotted histograms correspond to the contribution of individual sub-processes, i.e., q * +q * → γ * /Z → l + + l − , q * + g * → γ * /Z → l + + l − + q and q * +q * → γ * /Z + g → l + + l − + g, respectively. The shaded bands indicate the corresponding uncertainty ( 1 2 ≤ ζ ≤ 2) due to the hard scale variation with KMR UPDF in cross sections evaluation. Unlike the present report, most of the previous phenomenological works with the semi-KMR UPDF for the DY differential cross sections did not present the contribution of each sub-process in their final results and also, did not take into account the contribution of the third sub-process, assuming the possible double counting [50,51] between the first and the third sub-processes. However, according to our previous reports [25,53], we do not believe that there is any double-counting among the first and the third sub-processes. This point will be discussed in section III C.
(2): In the right panels of each figure, the results of different UPDF schemes applications, namely KMR, LO-MRW and NLO-MRW, in the differential cross sections are shown by the solid, dash and dotted-dash histograms, respectively, for the possible comparisons.

C. Discussions
First, we generally start by analyzing the calculated cross sections related to the medium and high center of mass energies for √ s = 1.8 and 7 TeV. Although the results show that the KMR UPDF describe reasonably the wide range of data of Tevatron and LHC, but for the two sets of differential cross section data which are in terms of p T and φ * η parameters, this is not the case. Indeed, in these two cases, the input NLO-MRW UPDF describe the data better than other schemes. To be sure about this conclusion, we try to include the newer data from the ATLAS collaboration at We should point out here that in our previous works, e.g., the reference [53], because of the possible fragmentation effects, the KMR and LO-MRW had a better agreement to the data.
The results of double and single differential cross sections dσ dM dy and dσ dM versus the invariant mass of the dilepton are compared to the experimental data at band of KMR, respectively. The experimental data also pass through the AOC band, see the reference [53]. A comparison between our results and the parton level Monte Carlo programs such as PYTHIA and SHERPA are also presented. According to these panels, although our results show a overestimate and underestimate in the low and high dilepton invariant mass region, the uncertainty band of our calculations covers the experimental data.
In the figures 5-8, the normalized differential cross-sections of DY at LHC as a function of the variable φ * η and different experimental conditions on the dilepton rapidity and invariant mass are presented. According to these figures (beside the figure 5), it is clear that in the small φ * η region, which corresponds to the back-to-back leptons, the contribution of the first and the second sub-processes are dominated and approximately in the same order. But in the figure 5, which is demonstrated for the small mass interval, all of the three sub-processes are in the same order for the small φ * η region. On the other hand for φ * η > 0.01 the three subprocesses are separated and as it is expected contribution of the first sub-process becomes enhanced and dσ 1 In the right panels of these figures, the same conclusion as above can be made about the effect of different UPDF schemes in which up to φ * η < 0.1 they behave the same, and for larger φ η as we discussed before, the NLO-MRW UPDF cross section calculations predict closer results to the corresponding data. The AOC and uncertainty bands approximately cover the ATLAS collaboration data. It is notable that the redefined form of normalization equation as xa(x, µ 2 ) µ 2 dk 2 t f (x, k 2 t , µ 2 ) without the factor 1/k 2 t does not lead to the collinear form of cross section after integrating over k 2 t . In the reference [51] the unpolarized DY in the pp collisions is investigated at the LHC energies by CCFM and semi-KMR within the reggeized quark formalism [51,52] to be sure about the gauge invariance of matrix elements. However, as we discussed in our previous work [53], the gauge invariance is guaranteed because of applying the small-x-approximation in our calculations. As we pointed out before, in the figures 1 (panel f) and 2 (panel b), our results are compared with those of references [51,52]. On the other hand, our results are compared with those of PYTHIA [79] (figure 1 (panels c-f)),

SHERPA [78] (figures 1 ( panels e-f) and 3 ( panels a,c,e)), FEWZ [80] (figure 7 (panels g-h))
and RESBOS [81] (figures 4 (panel e) and 7 (panel e)). It is observed that in the regions in which the higher order calculations are not important our results are similar to those of pQCD. However, in some parts, spatially for high p T and φ * η , the NLO-MRW UPDF scheme shows slightly different behavior with respect to the data and the pQCD methods.
In several papers such as the reference [77], the authors denote that the description of two observables, including the p T and φ * η distributions, are improved, if the higher order perturbative contributions are taken into account, which is in agreement with the cross check we performed. On the other hand as the scale of energy increases, for the LHC energies, the q-g sub-process has the largest contribution to the differential cross section in the most intervals of p T and φ * η , as it is expected. We also checked the interference effect between the γ * and Z in the cross sections and find out that the interference is ignorable in all regions.
Finally, we would like to point out that there is a new CMS measurement on the differential cross sections of the Z boson production in the P-P collisions [85]. In the figures 7 (in terms of p Z T ) and 8 (in terms of φ * η of dilepton) of this report, the CMS data are compared to the theoretical works presented in the reference [86][87][88][89], in which the UPDF (the so called transverse momentum dependent distribution functions (TMD)) are calculated, using the Parton Branching (PB) model. In this PB TMD model, the resummation to NLL accuracy, the fixed-order results at NLO, and the nonperturbative contributions are taken into account. The PB TMD results can predict the data well at low p Z T , but deviates from the measurements at high p Z T , because of missing contributions from Z+jets matrix element calculations. Furthermore, in the present work, we do not use the LHAPDF [90] or TMDlib [91] repositories, but we hope in our future reports, we can analyze the difference between the applications of present PDF and UPDF with those can be generated through LHAPDF and TMDlib repositories.

IV. CONCLUSIONS
We investigated the lepton pair production in the p-p and p-p collisions within the framework of k t -factorization approach. We used the transverse momentum dependent parton distribution functions of three different prescriptions, i.e., KMR, LO-MRW and NLO-MRW.
We calculated the matrix element square for the three different sub-processes among which the matrix element square for the q-q in the NLO level is rarely taken into account. We calculated several differential cross sections in terms of the dilepton invariant mass, transverse momentum and rapidity, as well as the angular correlation between produced leptons of the Drell-Yan process. In addition, we obtained the uncertainty band for the cross section distribution in the case of KMR by changing the scale factor as described in the section III. We considered the contribution of each sub-processes separately based on the off-shell and massless quarks. We found that although some of the results show that using the KMR framework, rather than LO-MRW and NLO-MRW schemes, represents more agreement with the experimental data, in the case of p T and φ * η probing the NLO-MRW gave better predictions. It is shown that the AOC and SOC constraints give similar results and our direct calculations of the off shell matrix elements and the method of integration for evaluation of the cross section give the same prediction as those of KaTie parton-level event generator.
Finally, in this work we consider the renormalization and factorization scales to be equal, i.e., µ R = µ F = ζM , in which, M is the invariant mass of produced dilepton and ζ can vary from 1/2 to 2, to estimate the scale uncertainties of our calculations. However as stated in the reference [85], one can vary each scale independently. Beside this it is possible to find the uncertainty, which come through the implementation of PDF through UPDF. But it should not be as large as the uncertainty effect due to the variation of renormalization and factorization scales. We hope to verify these effects in our future reports.  The matrix elements of three sub-processes can be presented as follows: where s = (k 1 + k 2 ) 2 and the electron and quark (fractional) electric charges are denoted by e and e q . Other notations are the same as the reference [50].

Appendix B:
The NLO splitting functions are defined as [21]: where i = 0 and 1 stand for the LO and the N LO, respectively. ∆ can be defined as [16]: and we have: panels. The contribution of q * +q * → γ * /Z → l + + l − , q * + g * → γ * /Z → l + + l − + q and q * +q * → γ * /Z + g → l + + l − + g sub-processes are presented by dash, dotted-dash and dotted histograms. In the right panels, the results corresponding to the KMR,          Uncertainty q q → l + l -          The differential cross-section of Drell-Yan lepton pair production at LHC as a function of the dilepton transverse momentum at E CM = 8TeV compared to the ATLAS data [33].The notation of all histograms is the same as in the figure 1.       as a function of φ * η at E CM = 8 TeV compared to the ATLAS data [33]. The otation of all histograms is the same as in the figure 1.  The normalized differential cross-section of Drell-Yan lepton pair production at LHC as a function of the dilepton rapidity at E CM = 7 TeV compared to the ATLAS data [84].
The notation of all histograms is the same as in the figure 1. LHC as a function of the dilepton rapidity at E CM = 7 TeV compared to the ATLAS data [84]. The notation of all histograms is the same as in the figure 1 .