Investigation of the Z boson production via hadron–hadron collisions in the k t and ( z , k t ) -factorization frameworks

In this paper, we study the Z boson production via the proton–proton (p–p) collisions within the k t and ( z , k t ) -factorization frameworks, using the Martin– Ryskin–Watt (MRW) unintegrated parton distribution functions (UPDFs) and the double unintegrated parton distribution functions (DUPDFs), respectively. For calculation of the differential cross section (DCS) within the k t -factorization ( k t is the partonic transverse momentum), the KATIE parton level event generator is used, while for the ( z , k t ) - factorization, the DCS is directly computed. Up to the tree level partonic next-to-leading order (NLO) are included, beside the inclusion of branching ratios, in our calculation. It should be noted that Martin, Ryskin and Watt are originally calculated the same process, i.e., the Z boson production, within the ( z , k t ) -factorization framework, while including only the lowest order tree level partonic sub-process. However, the present report extends their work from two perspectives. First, the additional sub-processes are included, second, for the processes up to the next-to-leading order (NLO), a direct calculation by considering the ﬁnal state leptons is imposed. Finally, we compare our results with the 13 T eV data of the ATLAS, LHCb, CMS collaborations, the corresponding collinear factorization predictions and the Modar-res, et al. reports. Our p–p DCS calculations show that the k t and ( z , k t ) -factorizations frameworks give relatively the same behavior in the central rapidity regions. While at the large rapidity regions, the ( z , k t ) -factorization, predicts the p–p DCS closer to the experimental data with respect to those of k t -factorization framework.

the complex nature of proton, where the partons are bound by the strong gluonic force inside the proton, this calculation is not an easy task. Hence, one cannot calculate the cross section, as simple as, that of the point like particles cross section. Fortunately, the factorization theorem allows us to calculate the cumbersome hadronic differential cross section (DCS), in terms of the elementary partonic cross section, convoluted by the familiar parton distribution functions (PDFs).
In the literature, there can be found different factorization frameworks for calculating the cross section, mostly based on the assumptions of parton behavior inside the proton, such as the collinear, k t , and (z, k t )-factorizations (k t is the partonic transverse momentum). For example, in the collinear factorization, it is assumed that the parton moves collinear to the proton and enters into the hard collision. Therefore in the collinear factorization theorem, it is supposed that the constituent partons can only have a fraction of the proton momentum, i.e., x, and these partons do not carry any transverse momentum dependency, i.e., k t . To describe the distribution of these partons inside the proton, one should use the parton distribution functions (PDFs), f a (x, μ 2 ), where μ 2 is the factorization (hard) scale. The scale dependency of the parton follows from the Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution equation [1][2][3], which describe the evolution of the collinear parton along the evolution ladder.
However, in the small x limit, the transverse momentum of parton (k t ), becomes comparable against the collinear component, i.e., x P (P is the proton momentum). Therefore, the evolution equation of collinear parton, should be generalized by the fact that, parton can also have transverse momentum. This means that the PDFs should have k t dependency, i.e., f a (x, k 2 t , μ 2 ). These k t dependent PDFs are called unintegrated parton distribution functions (UPDFs). Historically, due to the large role of the gluon at small x limit, UPDFs are defined for the gluonic content of proton. These UPDFs fol-low from the Balitsky-Fadin-Kuraev-Lipatov (BFKL) [4][5][6][7], and Catani-Ciafaloni-Fiorani-Marchesini (CCFM) [8][9][10][11] evolution equations. In fact the CCFM evolution equations by unifying the BFKL and DGLAP evolution equations are introduced, which allow to calculate the cross section in both small and large x regions. In order to achieve this goal CCFM uses angular ordering constraint along the evolution ladder. However, this model only works for gluon, and one cannot calculate UPDFs which also applies to quarks, which are necessary at large x limits. Therefore, other methods like Kimber-Martin-Ryskin (KMR) [12], Martin-Ryskin-Watt (MRW) [13], parton branching (PB) [14,15], etc are introduced to fill this gap. These UPDFs models allow to generate UPDFs in both relatively small and approximately large x regions, and hence one can apply the k t -factorization method in those limits.
The MRW UPDFs can be obtained by utilizing the DGLAP evolution equation. It assumes that the parton follows the DGLAP evolution equation to the last step of evolution ladder. But in the last step, parton becomes k t dependent by implementing the angular ordering on the gluon emission term. After that, the parton moves to the hard scattering without emitting any real emissions and becomes μ dependent. The MRW model is defined at both the LO and NLO levels. This model at the LO level is an extension of the KMR approach, wherein the angular ordering constraint is imposed on both the quark and gluon emission terms. In the Refs. [16][17][18][19][20][21][22][23], a detailed discussion, and investigation of these UPDFs are presented, and also for the applications in hadron (electron)-hadron collisions one can consult with the Refs. [23][24][25][26]. Beside the KMR and MRW methods for obtaining the UPDFs, another method is also recently introduced, which is called the parton branching (PB) approach, in which the UPDFs within this model can be obtained by iteratively solving the DGLAP evolution equation and collecting the transverse momenta of partons along the evolution ladder. One can see the applications of these UPDFs, for example, in the Refs. [23,27,28]. Among the different UPDFs discussed above, in this work, we only use the LO-MRW UPDFs, to perform the Z boson production DCS phenomenological study.
In addition to the factorization frameworks, which were discussed so far, another formalism is also introduced, that tries to generalize the k t -factorization by utilizing the full kinematics of the hard incoming and last step emitted partons. This framework is called the (z, k t )-factorization, and is introduced by Martin, Ryskin and Watt [29] (MRW DUPDFs). Previously, the DCS of different processes, such as, the electron-proton inclusive jet and diject productions [29,30], and the electroweak (W , Z ) bosons productions, in addition to the Higgs boson production [31], were calculated within this framework. In the (z, k t )-factorization approach, the last step emitted parton plays an important role in the DCS calculation. In contrast to the k t -factorization frame-work with the KMR and MRW UPDFs, where k 2 = −k 2 t , in this approach k 2 = −k 2 t /(1 − z), where z is the fractional momenta of hard parton, with respect to the parent parton. While in the k t -factorization approach, it is assumed that we are limited to the small x and z regions, in the (z, k t )factorization the full kinematics is considered. Therefore, the z dependency comes directly into the cross section calculation, and one should use the double UPDFs (DUPDFs), i.e., f a (x, z, k 2 t , μ 2 ), instead of the UPDFs. These newly introduced DUPDFs, as it will be discussed in this report, can simply be obtained via the MRW UPDFs. It should be mentioned that, because of this additional dependency on z, the calculation within this framework, becomes much more complicated with respect to other frameworks. However, despite these complications, it is interesting to test this framework in our DCS Z boson production phenomenological study, to gain more insight about this approach.
In this work, it is intended to directly calculate the Z boson production DCS for the first time within the (z, k t )factorization and compare its result to the k t -factorization approach, as well as, the collinear frameworks and the experimental data. The same process was previously calculated preliminary by Martin-Ryskin-Watt [31], in which they compare their results to the experimental data of C DF [32] and D0 [33]. However, their calculations suffer from two flaws. First, they do not consider the full kinematics of final state leptons, and calculate the cross section of Z boson production with the help of the branching ratio. Second, their calculations are only limited to the LO q + q → Z → l + l − sub-process. So, in this report, it is tried to fix these flaws by considering the complete final state leptons, including higher order sub-process, i.e., q + g → Z + q → l + l − + q, into the calculation. Additionally, in order to compare the above calculation with those of Ref. [22], we include the higher order sub-processes of this reference, with its corresponding parameters, in which the branching ratio factors are considered. Therefore, by doing this calculation, it is intended to perform a more detailed analysis of the k t and (z, k t )factorizations for the Z boson production compared to the original work of the Ref. [31]. Additionally, it should be mentioned that a similar analysis in the k t -factorization for the Z boson production is also performed in the Refs. [34,35], where the authors obtain good predictions of the data with the parton branching UPDFs approach.
The structure of this paper is as follows: in the Sect. 2, we first give a summary of the k t -factorization framework, which follows by the Sects. 2.1 and 2.2, to describe the MRW UPDFs and also numerical methods for calculating the pp DCS Z boson production, using the KATIE parton level event generator [36], respectively. In the Sect. 3, the (z, k t )factorization will be discussed, and furthermore in its three Sects. 3.1, 3.2, and 3.3, we explain the LO MRW DUPDFs, the partonic cross section within the (z, k t )-factorization approach, and also numerical methods within this framework, respectively. In the Sect. 4 our results and discussion are presented. The double counting issue will be described in the Sect. 5. Finally, the Sect. 6 is devoted to conclusions.

The k t -factorization framework
In the k t -factorization, one can write the p-p cross section as a convolution of the partonic cross section and the UPDFs. As it is mentioned before, in contrast to the collinear factorization framework, where it is assumed that the transverse momentum of parton has a negligible role in the above process, within the k t -factorization framework, the transverse momentum of hard parton has an important contribution. Therefore, one can generalize the collinear factorization formula, and write the p-p cross section in the k t -factorization as: where in the above formulaσ * ab , are the off-shell partonic cross sections, and f a(b) (x, k 2 t , μ 2 ) are the UPDFs. In this work, the MRW method is used for obtaining the UPDFs [13], and also in order to calculate the Eq. (1), we utilize the KATIE parton level event generator [36]. Henceforth, in what follows we discuss the MRW UPDFs, and also present a brief discussion about the KATIE parton level event generator.

The MRW UPDFs
The MRW UPDFs are introduced by Martin-Ryskin-Watt [13], which give the momentum distribution of partons inside the proton for both quarks and gluons. As it was pointed out in the introduction, in contrast to the BFKL and CCFM evolution equations, which are limited only to the gluon distributions, the MRW UPDFs do not have this restriction. This is due to this fact that the MRW UPDFs follow from the DGLAP evolution equations, and hence one can obtain the UPDFs for both quarks and gluons. In this approach, it is assumed that parton evolves up to the last evolution step collinearly along to the proton, obeying the DGLAP evolution equation, and only at the last step by emitting a real emission, the parton becomes k t dependent. Finally, parton evolves to the factorization scale μ without any real emission, using the Sudakov form factor. Therefore, one can write the MRW UPDFs as follows: where in the above equation, T a (k 2 t , μ 2 ) are the Sudakov form factors, f b ( x z , k 2 t ), are the momentum weighted the parton densities at the LO level, which can be either for quark (antiquark) or gluon, respectively, and δ ab (z max − z) is constraint on the gluon emission, when a = b. The Sudakov form factor has the following forms: In order to obtain the cutoff on z and ξ (to constraint the gluon emission in the Eqs. 2, 3), one can either use the strong or angular ordering in the last evolution step. The strong ordering cutoff (SOC) is in fact an approximation of the angular ordering cutoff (AOC) in the large z limit, i.e., z → 1 [37]. Therefore, if one is interested in the effect of small x limit, it is better to use the AOC. The AOC cutoff is written as: where q = k t /(1 − z), and in the last step of the evolution equation, the above cutoff can be written as: If one considers the SOC, i.e., the z → 1 limit, the above cutoff can be written as: It can be seen that the SOC cutoff is tougher and limits the parton's transverse momentum to the k t < μ regions, while the AOC is smoother and the parton transverse momentum can be larger than the factorization scale. As it can be seen from the Refs. [38], the Kimber-Martin-Ryskin (KMR) approach was first introduced with the SOC. However, as it is stated in the Refs. [39,40], the SOC cutoff is replaced by the more correct version AOC. Later, the Martin-Ryskin-Watt (MRW) approach only used AOC [13]. Therefore, due to what was discussed above, we also adopt AOC in this work. However, in general, the angular ordering form of this cutoff is more preferred in the literature on the KMR and MRW approach, and hence in this work we adopt this cutoff for the MRW UPDFs which is, z max = μ (μ+k t ) and ξ max = μ (μ+κ t ) . It should also be mentioned that in the Ref. [13], it is stated that the above integral version of MRW UPDFs, can also be expressed as the following "differential form": However, it is shown in the Refs. [41][42][43] that, the integral and differential versions of this UPDFs, do not give the same UPDFs. In fact, using the differential version is problematic, and in some parton transverse momentum regions, it gives us the negative and discontinuous UPDFs. Therefore, as it is demonstrated in the Ref. [41], in order to obtain the same UPDFs as the integral version, in the differential form, the cutoff dependent PDFs and an additional term should be used. Another aspect of the Ref. [43] is the question of the KMR UPDFs validity, and it claims the wrong application of the modified DGLAP evolution equation. In this work it is pointed out that the virtual term of the modified DGLAP evolution equation lacks a 'z' factor, i.e.: and it leads to a wrong form of Sudakov form factors. However, by considering the results of Refs. [39,40], it is clear that KMR uses correct Sudakov form factors in their work. In other words, the modified DGLAP evolution equation is implemented in a condensed form, but in a confusing manner. Therefore, it is evident that the reasoning of Ref. [43] on the validity of the KMR approach is not justified. Actually, as it is stated in the Ref. [13] the difference between the KMR and MRW approach is on the different imposition of AOC. In fact in the KMR, the AOC is imposed on both the quark and gluon emission terms, while in the MRW, the AOC is only imposed on the gluon emission term. In the Ref. [43], also one of the challenges in the non-equivalency between the differential (DKMR/DMRW) and the integral form of KMR/MRW (IKMR/IMRW) at k t > μ is solved. Generally, we should note that, DKMR/DMRW can be problematic, and it can give us the negative distribution especially at large x and k 2 t regions. However, by adding the additional term that is obtained, both in our work [41] and also with different approach in the Ref. [43], one can obtain acceptable UPDFs. It should also be mentioned that in the work of Ref. [43], it is claimed that with this new formula one can obtain the same UPDFs from both of the DKMR/DMRW and IKMR/IMRW UPDFs. But, in our work [41], we showed that one needs also the cutoff dependent PDFs in DKMR/DMRW, as proposed in the Ref. [42]. The important point here is this fact that cutoff dependent PDFs only become important at large x regions, when the cutoff becomes much important. However, it should be noted that we use the integral version of the KMR/MRW formalism, so our work is consistent with the findings of the Refs. [41][42][43].
It should be mentioned that the Eqs. 2, 7 are only applicable at k t ≥ μ 0 , with μ 0 ∼ 1 GeV, since the input PDFs, f a (x, k 2 t ), are not defined at the scales less than μ 0 . As a result, to obtain the MRW UPDFs at k t < μ 0 , one can conveniently utilizes the normalization constraint, i.e.,: and obtain the MRW UPDFs, that satisfy the normalization condition. On the other hand, for k t < μ 0 , it can be shown that the following distribution can satisfies the normalization condition: 1 However, it is necessary to point out that, even by using the above distribution, which is also introduced in the Ref. [13], the MRW UPDFs still does not satisfy the normalization condition. Because, by looking at the derivation of the Ref. [13], one notice that the Eq. 10 can be obtained by using the differential version of the MRW instead of the integral one of these UPDFs. Therefore one task, that would be of interest, is to obtain the appropriate MRW UPDFs for k t < μ 0 , which can fully satisfy the normalization condition, and we leave it off for our future works. Eventually, the LO MRW UPDFs for the (anti) quarks and gluon and their corresponding Sudakov form factors can be written as follows: with the Sudakov form factors as: where in the above equations q is denoted for u, u, d, d, . . ., and also in the Eq. (14), n F is the number of active quark and anti-quark flavors and ξ min = 1 − ξ max . Because we compare our predictions with the results of Ref. [22], in which their calculations are performed with the KMR UPDFs, it should be mentioned that the KMR and MRW UPDFs at the LO level do not have a significant difference in their definitions. But since, in the MRW UPDFs the angular ordering cutoff is only imposed on the soft gluon emission terms, while in the KMR UPDFs this cutoff is imposed on both gluon and quark (anti-quark) emission terms, where this cutoff ordering, actually leads to the some difference between these two approaches only in the large p ll t and x values, where this cutoff becomes important.

Numerical methods
In the present report for calculating the cross section within the k t -factorization framework, the KATIE parton level event generator [36] is utilized. With this parton level event generator, one can obtain the partonic cross section, as well as, the DCS with respect to different variables at the tree level for both the collinear and k t -factorization frameworks. There are two different approaches for providing UPDFs as inputs to this parton level event generator. One should either use TMDlib [44], (which for example generates the differential form of MRW UPDFs), or in the case where the desired UPDFs are not available in this library, the UPDFs grid files should directly be provided. Because the integral version of MRW UPDFs is not available in the TMDlib library, for present calculation, the MRW UPDFs grid files are directly generated through our program file, which is provided publicly at the Ref. [23]. For the input PDFs of the MRW UPDFs in the Eq. (2), the MMHT2014lo68cl [45] set of the LHAPDF [46] library is used.
In the KATIE parton level event generator, we can calculate the Z boson production cross section by considering the five sub-processes, i.e., q + q → Z → l + l − , q + g → Z + q → l + l − + q, g + g → Z + q + q, g + q → Z + g + q and q + q → Z + g + g. However, the first 2 sub-processes have the main contribution with respect to the higher order ones, i.e., q + q → Z → l + l − and q + g → Z + q → l + l − + q. When calculating these sub-processes, three points should be considered. First, in order to avoid soft gluon divergence, we set a cutoff on the transverse momentum of final state gluon to be larger than the factorization scale (i.e., the final gluon to be hard). Second, the KMR and MRW UPDFs have limited transverse momentum in their definition, on the other hand, there is the Sudakov form factor, that have a suppressive effect and can make the cross section well-behaved and finite. Additionally, in our previous work [23], it was shown that there is no divergence indeed occurring in the k t -factorization formalism. This was obtained by calculating the differential cross section with respect to the rapidity of final state quark in the q + g → γ γ γ + q sub-process. It was demonstrated that the differential cross section is well-behaved in the large rapidity limits. Third, it is also claimed that there is a double counting between 2 → 2 and 2 → 3 sub-processes. However, as it is discussed in the Sect. 5, this topic is controversial, and even if the double counting is considered, it does not have huge impact in our results.
We perform these calculations for the first three quark flavors, i.e., up, down, strange, and their anti-quarks. μ f = μ r = p ll 2 t + m ll 2 as the factorization and renormalization scales is chosen, in which p ll t and m ll are the transverse momentum and the invariant mass of the output dilepton, respectively. Also, we consider the experimental cuts of ATLAS [47], LHCb [48], and CMS [28] experiments in our calculation with this generator, which are expressed in detail in the Sect. 4.

The (z, k t )-factorization framework
In contrast to the k t -factorization, where the full kinematics of hard parton is not considered, in the (z, k t )-factorization the full kinematics is taken into account. Within the (z, k t )factorization, there is an additional component, i.e., light cone minus term, where it can be fixed by the constraint of on-shellness of the last step emitted parton. As a result of this, in the (z, k t )-factorization we have z, in addition to the k t dependency. Therefore, considering the complete kinematics of the parton at the last evolution step within the (z, k t )factorization, in contrast to the k t -factorization approach, one can expect better results in the forward rapidity limits, wherein the last step parton becomes important [29,30]. This additional degree of freedom, also makes the incoming hard parton virtuality to become k 2 = −k 2 t /(1 − z) within this framework, while in the case of k t -factorization the above virtuality is k 2 = −k 2 t . Another important effect of this additional degree of freedom, i.e., the light-cone minus term in the (z, k t )-factorization, is this fact that within this formalism the UPDFs become z dependent, and hence one needs to modify hadronic factorization formula for calculating the cross section, compared to the k t -factorization approach. Therefore, by generalizing the k t -factorization framework, one can write the general p-p cross section formula in the (z, k t )factorization as: whereσ * ab are the off-shell partonic cross sections and f a (x, z, k 2 t , μ 2 ) are called double UPDFs (DUPDFs). As it can be seen in the above formula, the DUPDFs and also the partonic cross section, in addition to dependency on k 2 t and μ 2 , they have a dependency on z. Therefore, in order to calculate the hadronic p-p cross section, one needs to devise methods for calculating both DUPDFs and also the partonic cross section, which is in accordance with this new kinematics.

The LO MRW DUPDFs
The DUPDFs are first introduced as a generalization of the MRW UPDFs model by Martin, Ryskin, and Watt [29]. In order to obtain these DUPDFs, one needs to simply ignore the integral over z of the MRW UPDFs approach. As a result of this, we can write the DUPDFs similar to the MRW approach, i.e., the DMRW model, for both quark, (anti-quark) and gluon, respectively, as follows: where in the above equations q stand for u, u, d, d, . . . and also z max , ξ max , n f , and ξ min are defined below the Eq. (3) and the Eq. (14), respectively. Additionally, one should note that in order to satisfy unitarity, these DUPDFs like UPDFs should also satisfy the normalization condition, which has this modified form:

The partonic cross section within the DMRW model
Another ingredient of the cross section formula, i.e., the section 15, within the DMRW framework is the partonic cross section. In the Fig. 1, DMRW approach with its corresponding kinematics and also its relation with the partonic cross section is schematically shown. As it is shown in the Fig. 1, the two incoming hard partons separately emit a parton with the momentum k i = l i − k i , where i = 1, 2, before the hard collision. After their last step emissions, these partons with the momentum k i collide with each other. Due to the nature of the collision and for simplicity in performing the calculation, we use the Sudakov decomposition of momenta of two incoming partons [31,49], i.e.: where i, j = 1, 2 or 2, 1. Additionally, we also neglect the hadron masses in the center of mass frame of the colliding hadrons. Hence, the squared center of mass energy can be written as s = (P 1 + P 2 ) 2 2P 1 ·P 2 . So one finds the following relations: In the above equations, the definition of a = (a + , a − , a ⊥ ), with a + = (a 0 + a 3 ) and a − = (a 0 − a 3 ) for the light-cone variables is used. In the Fig. 1, the penultimate parton in the evolution ladder has momentum l i = (x i /z i )P i and hence using the Eq. (21), the momentum of the last step emitted parton along the evolution ladder is obtained, i.e.,: Fig. 1 Illustration of the DMRW model in the hadron-hadron collision. In the left panel of this figure, the transverse momentum of each incoming parton into the subprocess is generated by a single parton emission in the last evolution step. While, in the right panel of this figure, the last evolution step is factorized into the DUPDFs, i.e.,

Fig. 2
The left and middle panels Feynman diagrams correspond to the q + g −→ Z + q −→ l + + l − + q processes and the right one corresponds to the q + q −→ Z −→ l + + l − process The β parameter in the Eq. (21) is determined by the on-shell condition for the emitted parton, i.e., k 2 i = 0: where r i ≡ k 2 it /s. Using the above kinematics in the DMRW model, one can show that the virtuality of two incoming partons is k 2 i = −k 2 it /(1 − z i ). This is in contrast to the k tfactorization, where z = 0 and β = 0, with the virtuality k 2 = −k 2 t . For calculating the cross section within the DMRW approach, the following two sub-processes are included: We first start with the cross section calculation of the Sect. 15, with the sub-process of the Eq. (25), which is also represented schematically in the right panel of the Fig. 1. Before proceeding the following discussion, it should be mentioned that the above sub-processes are considered directly into our calculation, by taking into account the full kinematics (Fig. 2). The partonic cross section in this figure, σ q 1 q 2 , can be written as: where for the sub-process of the Eq. (25), q 1 (q 2 ) is q(q), the F q 1 q 2 = 2x 1 x 2 s is the flux factor, q 1 q 2 denotes the phase space integrals and M q 1 q 2 is the amplitude. The amplitude of the sub-process in the Eq. (25), can also be calculated with the help of Feynman rules and written as follows for this sub-process: where m Z and Z are the mass and the full decay width of Z boson, θ w is the weak mixing angle,ŝ is the partonic center of mass energy and C V and C A are the vector and axial constants.
In order to calculate the cross section for the sub-process of the Eq. (25), we should calculate: where φ 1 and φ 2 are the azimuthal angles of the incoming partons. For calculating the above integrals, one needs to simplify the above Dirac delta functions, corresponding to the 2 → 2 conservation of energy and momenta: where in the above equations y 3 and y 4 are the rapidities of the final state leptons. So by imposing the Eq. (30) one can remove the k 4t integral in the Eq. (29), and also using the Eqs. (31) and (32), the above mentioned integrals with respect to x 1 and x 2 can be simplified. Due to the complicated relations for x 1 and x 2 , i.e., Eqs. (31) and (32), they will be calculated with the help of the Mathematica software. Accordingly, the following relation for x 1 can be obtained: where A 1 in the above equation is return as: × e y 4 k 3t + e y 3 k 4t e 2y 3 +y 4 k 3t + e y 3 +2y 4 k 4t . (34) In the Eq. (33), only x − 1 gives us the physical result, and x + 1 does not play any role in our calculation. Hence in what follows, we write x − 1 as x 1 for simplification. Having x 1 obviously the following equation can be obtained for x 2 : After, simplifying the Eq. (29), we can find the following relation for the cross section of the 2 → 2 sub-process: Then, similar to the above 2 → 2 sub-process, one can also write the cross section for the 2 → 3 sub-process, i.e., Eq. (26), as follows: where φ 1 and φ 2 are the azimuthal angles of the incoming partons, while φ 3 and φ 4 are the azimuthal angles of the first and second outgoing particles. The amplitude of the subprocess in the Eq. (26) can also be calculated by the Feynman rules as: The integrals in the Eq. (37) can also be simplified using the Dirac delta functions corresponding to the following conservation of energy and momenta relations: After performing simplifications in a similar manner as what we did for the Eq. (29), with the help of the Mathematica software, the relations for x 1 and x 2 can be obtained. Due to the lengthy and complicated relations for x 1 and x 2 , their full expressions are not expressed here. Performing these calculations similar to the 2 → 2 sub-process, the Eq. (37) becomes as follows: Additionally, to compare our calculations with those of the Ref. [22], we calculate the cross section of higher order sub-processes, i.e, g + g → Z + q + q, g + q → Z + g + q and q + q → Z + g + g, using the branching ratio method. The Feynman diagrams of these sub-processes are presented in the Fig. 3. We also adopt f (Z 0 → e + + e − ) = 0.0336 [50], as the branching ratio in our calculation. This choice is mostly due to the increase of numerical speed and the reduction of integrals. Also, it should be stressed that these sub-processes, as it will be shown in the Sect. 4, have negligible roles in the final results. Because these sub-processes have the kinematics of 2 → 3, therefore one can use the Eq. (42), but by replacing with the parameters involved in the three aforementioned sub-processes, including their matrix elements, as shown in the Appendix 1. Additionally, similar arguments for double counting and divergency still hold in the (z, k t )-factorization as in the k t -factorization. Finally, it should be stressed that to avoid soft gluon divergence, we set a cutoff on transverse momentum of final state gluon to be larger than the factorization scale like what we performed for the k t -factorization. On the other hand, there is also another divergence when k 2 t → 0, i.e., final state jets become collinear to the incoming hadron. But, in accordance to the Ref. [23], this divergence does not happen in the k t -factorization. Because, our UPDFs and DUPDFs are convoluted to the partonic cross section, and they are set to zero when their transverse momenta are less than 0.03, i.e., f a (x, k t < 0.03, μ) = 0 and f a (x, z, k t < 0.03, μ) = 0. It is also shown in the figure 11 of the Ref. [23], that with this aforementioned constraint on the minimum transverse momentum of the incoming parton, results are convergent within the k t -factorization cross section calculation. After this lengthy discussion about the DMRW model, in the following section, we discuss our numerical methods for calculating the above integrals and also other parameters involve in our work.

Numerical methods within the DMRW model
Generally, the cross section calculation within the DMRW model, because of the additional integrals, and also its complicated kinematics is more difficult with respect to the MRW framework, wherein the full kinematics of hard parton is not considered. The main component of these two equations, i.e., the Eqs. (36) and (42), is the matrix elements squared of Eqs. (28) and (38), corresponding to the q + q → Z → l + + l − and q + g → Z + q → l + + l − + q sub-processes. The aforementioned matrix element squared is calculated, with the help of Feynart [51] and FeynCalc [52] Mathematica packages. It should also be noted that in order to be consistent with the method of Refs. [29,31], in calculating the DMRW approach amplitude, the momenta of incoming partons are only kept in their full kinematics in the denominator of the scattering amplitude, while in its numerator the collinear kinematics is kept.
To compute the cross sections, the multidimensional integrals of the Eqs. (36) and (42) are performed with the VEGAS Monte Carlo method. Additionally, similar to the MRW case, for the factorization and renormalization scales, we choose μ f = μ r = p ll 2 t + m ll 2 , where p ll t and m ll are the transverse momentum and the invariant mass of the output dilepton, respectively. Also, the non-logarithmic loop corrections to the sub-process of the Eq. (25) are taken into account by including the K -factor in our calculation, as: where C F = 4/3, and the particular choice μ 2 = p ll 4/3 t M 2/3 Z is proposed to eliminate the sub-leading logarithmic terms [31]. Finally, it should be mentioned that the MRW DUPDFs (DMRW) can simply be obtained similar to the MRW UPDFs, but here we do not perform the integration on z parameter. Also, to speed up our calculation only the first three quark flavors are taken into account, i.e., up, down, and strange.

Results and discussion
One of the useful tools to study the hadronic interactions in the high energies is the Z boson production. The intermediate Z bosons are mostly produced through qq annihilation decaying into the lepton pairs, i.e., the sub-process of the Eq. (25), and also from the Compton scattering of subprocess in the Eq. (26). In this section, it is intended to discuss the results of Z boson production within the MRW and DMRW factorizations by comparing these predictions with the 13 T eV data of the ATLAS [47], CMS [28], and LHCb [48], and also their corresponding collinear counterparts. As it was noted before, a similar calculation in DMRW model is performed originally by Watt, Martin, and Ryskin [31]. However, as it was pointed out in the introduction, their calculations are simplified from two sides. First of all, in their report, only one 2 → 2 sub-process is considered, i.e., the Eq. (25). Second, they do not fully consider the leptonic final state, and calculate their results with the help of the branching ratio of Z boson to the corresponding leptonic final state. While, our calculations intended to performed, by fixing these two flaws of the Ref. [31]. Additionally, to compare our results with the work of Ref. [22] (calculated by considering up to NNLO Fig. 3 The 8 left panel Feynman diagrams correspond to the q + q → Z + g + g, the 8 right panel Feynman diagrams correspond to the g + g → Z + q + q and the 6 bottom panel Feynman diagrams correspond to the g + q → Z + g + q sub-process sub-processes using the branching ratio), and also see the role of the higher order sub-processes, some of these higher order sub-processes corresponding to the aforementioned work is calculated with the help of the branching ratio [50]. Before presenting our results, it should be mentioned that in all of the figures, the MRW and DMRW predictions with the MRW UPDFs and DUPDFs, are all denoted by MRW and DMRW, respectively. Also, all collinear results presented in our work are extracted from the corresponding experimental papers and we only did use the result of MadGraph5-aMC@NLO, Pythia8 and ResBos event generators from mentioned references, respectively, [28,47,48]. Therefore, if one needs a detailed discussion of the calculation method, should be consulted with the reference of each experiment mentioned in this work. Finally, as it will be presented in our results, as one should expect, generally, the collinear factorization shows a better behavior with respect to the MRW and DMRW frameworks. However, as it is discussed in the Ref. [31], although the result of the MRW and DMRW frameworks are not as good as those of collinear factorization, but calculating cross section with the MRW UPDFs and DUPDFs is more simplistic, since they reduce the numbers of Feynman diagrams and the computational time with respect to that of collinear factorization.
In the Fig. 4 our predicted 1/σ dσ/d P ll t and 1/σ dσ/dφ * η , is compared with those given by 13 TeV ATLAS data. In this experiment, the pseudorapidity of each final state electron is limited to the |η| < 2.47, excluding 1.37 < |η| < 1.52. The invariant mass of the final state leptons is limited between 66 GeV < M e − e + < 116 GeV. Additionally, the transverse momenta of each final state electron are larger than 27 GeV. In the left and right panels of this figure, the above DCS with respect to P ll t , and φ * η in the k t and (z, k t )-factorization frameworks, i.e., MRW and DMRW approaches, are compared to the ATLAS data,as well as, the Pythia8 results [47]. It can be seen from this figure, that both of the MRW and DMRW results, show relatively the same behavior. This is consistent with this fact that they mostly become different from each other, at large rapidity regions, and also large z limit, where one can expect that the DMRW behaves better than the MRW. Despite this similarity, at small dilepton transverse momenta, the DMRW tends to be larger than the MRW. This behavior is actually due to the effect of the K -factor, i.e., the Eq. (43), which leads to this conclusssion that the DMRW become larger compared with the MRW in this region, where K -factor is omitted. Also it should be mentioned that, both of these frameworks tend to overestimate the data at large transverse momenta. However, it is obvious from this figure that the DMRW shows much better prediction at medium range of dilepton transverse momenta relative to the MRW. The reason for this overestimate, at large transverse momenta of lepton pair, is mostly due to the MRW and DMRW models. For example, if instead of the LO MRW, one uses the NLO-MRW or PB UPDFs models, a better result at the large Fig. 4 The Left (right) panel shows the comparison of MRW and DMRW results using the MRW UPDFs and DUPDFs with the 13 T eV ATLAS experimental data [47], and also the Pythia8 [47] results for the DCS with respect to p ll transverse momenta is obtained [23]. This issue is discussed in the k t -factorization framework, for example by comparing the results of the various DCS, using the LO MRW model with those of the NLO MRW or PB UPDFs, see the Refs. [23][24][25]. As it was explained in these references, the reason for this overestimation of LO MRW model with respect to those of NLO MRW and PB UPDFs is due to this fact that at k t > μ, the corresponding the UPDFs have a larger contribution compared with other UPDFs such as the NLO-MRW or PB UPDFs. For example in the case of the NLO-MRW, one has k 2 = k 2 t 1 − z < μ 2 virtuality ordering constraint, i.e., (μ 2 − k 2 ). This constraint within this model becomes 0 at k t > μ. Hence, one expects that the NLO MRW, as it is shown in the aforementioned references, to perform much better in those regions.
In the right panel of the Fig. 4, the MRW and DMRW DCS with respect to the φ * η are plotted. This variable is defined as: where φ acop = π −| φ|, and cos(θ * η ) = tanh[(η − −η + )/2]. Here φ is the azimuthal angle (in radians) between the two leptons, and η −(+) is the pseudorapidity of the negatively (positively) charged leptons. The parameter φ * η is more of interest to the experimentalists compared to the dilepton transverse momentum. Mostly, because this parameter dependens only on measuring the direction of final state leptons, rather than their momenta. By observing the right panel of this figure, in contrast to dσ/d P ll t , where the DMRW tends to be larger than the MRW, in the small dilepton transverse momentum, here in the range of small and medium φ * η , the exact opposite behavior occurs, and the DMRW underesti-mates the data relative to the MRW. Although, the DMRW gives better description of the data in the range of large φ * η with respect to the MRW. Additionally, as it was pointed out before, the Pythia8 results cover the data well in all regions.
In the Figs. 5, 6 and 7 the MRW and the DMRW for the 13 TeV LHCb data [48] are presented. Additionally, in the Fig. 5, a comparison with the ResBos [48] results is presented. The experimental data belong to the Z → μ + μ − process, where each of the muons is within the pseudorapidity interval of 2.0 < η < 4.5 and has the transverse momentum p T > 20 GeV. Additionally, the pair of produced muons invariant mass is required to be 60 GeV < M μμ < 120 GeV.
In the top left and right panels of the Fig. 5, a comparison between the DCS of KMR, MRW and DMRW, with respect to the P ll t and φ * η is shown, respectively. It is observed that similar to the ATLAS data in the Fig. 4, the DMRW gives much better description of data relative to the MRW in the range of small to medium transverse momenta. While the DMRW in the range of small φ * tends to underestimate the data, which is in contrast to the MRW, that can cover the data well within this limit. Additionally, by comparing the KMR prediction of the Ref. [22] (calculated with the branching ratio method), with the MRW, which is calculated with the KATIE event generator, one can notice that these two UPDFs, i.e., the KMR and MRW, are only similar to each other in the intermediate range of transverse momentum, and they become different from each other at the small and large transverse momenta. One of the main reason for these differences between the KMR and the MRW results can be due to this fact that, the different implementation of angular ordering cutoffs is utilized in this two models. However, since their calculation methods are different, other parameters may also be involved. Also, one should note that, the results of Ref.  [22], the MRW and the DMRW results with the 13 TeV LHCb experimental data [48], and those of ResBos [48] results with respect to p ll t (φ * η ). In the left bottom panel, the same comparison for DCS with respect to y Z is shown. The right bottom panel of this figure shows the relative contribution of the higher order sub-processes (the sum of g+g → Z +q +q, g+q → Z +g+q and q +q → Z +g+g), denoted by σ 2 , and lower order sub-processes (the sum of q + q → Z → l + l − and q + g → Z + q → l + l − + q), denoted by σ 1 , for the MRW and the DMRW [22] is limited to p ll t < 50 GeV, and the comparison with the data in all regions is not reported. The last thing to note in this figure is that, the ResBos results, like what was observed in the Pythia8 of the ATLAS experiment, give better description of the data compared to both of the MRW and the DMRW.
In the bottom left panel of Fig. 5, one observes the variation of KMR, MRW and DMRW approaches with respect to the y Z . It is evident from this figure that, their results are close to each other in the Z boson rapidity region of y Z < 4, while they become separate from each other in the y Z > 4, wherein the MRW fails to describe the data well in that region. This issue can be attributed to the problematic behavior of the MRW in defining the last step emitted parton, which is discussed at the beginning of Sect. 3. In fact, the last step emitted parton can actually play an important role in the forward region, hence one can expect that the MRW to be worse compared to the DMRW. Therefore it is essential in this region to treat correctly the last step emitted parton. Although it should also be stressed that the difference between these two frameworks essentially is not significant. One final note about this figure is that, interestingly the KMR UPDFs of the Ref. [22], can cover the data well in all y Z regions, including the forward ones. This may be due to the different method of calculation of Ref. [22], which are done with the branching ratio, while in this work our the MRW is performed with the KATIE parton level event generator. Additionally, despite the fact that the ResBos result, can cover the data well within all rapidity regions, but it tends to overestimate the data at large y Z limit.
In the bottom right panel of the Fig. 5, a comparison between the contribution of the higher order sub-processes (the sum of g + g → Z + q + q, g + q → Z + g + q and q + q → Z + g + g), denoted by σ 2 , and lower order sub-processes (the sum of q + q → Z → l + l − and q+g → Z +q → l + l − +q), denoted by σ 1 , for the MRW and the DMRW are compared. As it was discussed before, in the Ref. [22], the aforementioned higher order sub-processes are included. Therefore, their contributions in the dσ /dp ll T DCS are calculated. As it is obvious from this figure, the role of higher order sub-processes, i.e., those calculated in the Ref. [22], is negligibly relative to the lower order sub-processes. Therefore one can safely ignore their contributions into our calculation.
In the Fig. 6, it can be observed the double DCS with respect to p ll t in various rapidity regions of the produced Z boson. Similar to our previous results for the cross section with respect to p ll t , it can be obtained relatively the same behavior in all of the regions except where 4 < y Z < 4.5. As it is obvious from the left bottom panel of the Fig. 5, i.e., dσ/dy Z , it can be expected similar worse behavior for the MRW compared to the DMRW within the large rapidity regions. In fact as we move toward large rapidity regions, the DMRW becomes much better relative to the MRW, especially in small and large dilepton transverse momentum regions. Also, while at small rapidity region, the DMRW treats better especially at the small dilepton transverse momentum, the MRW tends to underestimate the data.
In the Fig. 7, the double DCS predictions of the MRW and the DMRW with respect to φ * η in the various y Z regions are shown. Generally, the results are in accordance with our expectations, and the DMRW underestimates the data at small φ * η , while its predictions cover the data at large φ * η . Also, in the most φ * η and y Z regions, the MRW can cover the data well. However, as it was discussed before, at large rapidity region, i.e., 4 < y Z < 4.5, the DMRW results behave better, and its prediction covers the data relatively well in all φ * η regions, while the MRW underestimates the data in all φ * η regions, even for those of small φ * η limit. In the Fig. 8, the results of MadGraph5-aMC@NLO [28], the KMR, the MRW and the DMRW for the 13 TeV, CMS data are shown. These data belong to the Z → l + l − process, where each of the leptons is within the pseudorapidity interval of |η| < 2.4. These produced leptons are required to have the transverse momentum p T > 25 GeV and the invariant mass of the pair of produced leptons is required to be |m ll − 91.1876 GeV| < 15 GeV.
Also, in the Fig. 8, the DCS with respect to P ll t , φ * η , and y Z are shown. Similar to our results of the MRW and the DMRW for the 13 TeV ATLAS and LHCb data, it can also be observed relatively the same behavior for the 13 TeV CMS data. It should be mentioned that because in this experiment the forward regions, i.e., 2.5 ≤ y z ≤ 4.5, do not play any role in calculation, therefore it is expected that, the MRW and DMRW become similar to each other for the ATLAS experiment. Comparing the predictions of KMR UPDFs of Ref. [22] with our results, also show that this model can cover the data of dσ/dy Z , and also dσ /d P ll t . It should be noted that again the result of the dσ/d P ll t is limited to P ll t < 50 GeV region, and therefore a full comparison with our work in all dilepton transverse momentum regions is not possible. Finally, it should be mentioned that the collinear results of MadGraph5-aMC@NLO can predict the data well in all of the channels in the Fig. 8.

The glance of the double counting
Before presenting conclusions, it is useful to discuss an important topic about the double counting between the lower and higher order sub-processes in the k t -factorization. Some authors argue that there is the double counting, for example, between q + q → Z → l + + l − and q + g → Z + q → l + + l − + q sub-processes, and if one address the double counting, the predictions will be improved at the large dilepton transverse momentum. However, as stated in the Refs. [53,54], there is no double counting between the 2 → 2 and 2 → 3 sub-processes. The KMR and MRW UPDFs in the k t and (z, k t )-factorizations should correspond to the PDFs in the collinear case. All the splittings and real emissions of partons, including the last step emission, are indeed factorized in the UPDFs. Additionally, any changes to the UPDFs have a problematic effect on the normalization condition that relates UPDFs to PDFs. Another argument that contradicts the idea of having double counting is also mentioned on page 536, paragraph 4 of Ref. [54]. However, we are also interested in checking how our results would be affected if we assume that there is double counting between the aforementioned sub-processes. In order to resolve the double counting issue, we utilize the method proposed in Ref. [55]. In accordance with the approach of this article, we limit the transverse momenta of incoming partons in the q + q → Z → l + + l − sub-process to the regions less than the factorization scale. Additionally, in the sub-process q + g → Z + q → l + + l − + q, we impose a cut on the final quark transverse momentum such that to be larger than the incoming partons transverse momentum. According to the Ref. [55], the first constraint excludes the possible extra hard emissions from the UPDFs, and the second constraint also limits the final state jet to contribute from the matrix element. We calculate our k t -factorization results in accor-dance with this approach, which is shown in the Fig. 9. As can be seen from the Fig. 9, even treating the double counting with this approach, the significant effect dose not appear in our results. It should be noted that the UPDFs alone have the biggest role in our predictions. Moreover, considering the large contribution of q + q → Z → l + + l − process(see the Fig. 9), we can observe that even this simple sub-process overestimates experimental data at large dilepton transverse momenta. Therefore, it is clear that even if it were assumed that there is a double counting, treating it, would not improve our results at large dilepton transverse momentum.  [22], MRW and DMRW models with the 13 TeV CMS experimental data [28], and also MadGraph5-aMC@NLO [28] results for the DCS with respect to p ll t (φ * η ). In the bottom panel, the same comparison for the DCS with respect to y Z is shown

Conclusions
In this work, the Z boson production cross section calculated within the k t and (z, k t )-factorizations frameworks using the MRW UPDFs and DUPDFs by considering q + q → Z → l + + l − and q + g → Z + q → l + + l − + q sub-processes. Furthermore, this study examines the role of some higher order sub-processes that are computed in the Ref. [22] by employing the branching ratio method and compares our predictions with the outcomes of this reference. Our results are compared with various channels of 13 T eV data of ATLAS [47], CMS [28], and LHCb [48], and also their corresponding collinear results. This calculation is performed for the first time within the (z, k t )-factorization by including one higher order sub-process and also keeping the kinematics of the final state dilepton into our calculation. It is showed that both of the MRW and the DMRW behave relatively similar to each other. However, we observed that the MRW UPDFs in the rapidity region of 4 < y Z < 4.5 fails to describe the experimental data. This is in contrast to the DMRW result, which is successfully cover this limit. The reason for this behavior is that the MRW, in contrast to the DMRW does not treat the last step emitted parton along the evolution ladder as a real detectable particle, and hence in this forward region, we can expect this framework to fail. Additionally, It is also observed that in the small and medium transverse momentum of final state dilepton, the DMRW shows better behavior with respect to its counterpart. However, it was showed that at small φ * η , the DMRW underestimate the experimental data, although the MRW is pretty successful in those regions. Also, a comparison with the results of the Ref. [22] shows that the KMR UPDFs describe the data of the LHCb and CMS data much better compared to the MRW and the DMRW models due to the additional constraint on the quark emission terms Fig. 9 Comparison between different contributions of the q + q → Z → l + + l − , and q + g → Z + q → l + + l − + q sub-processes, calculated based on the method proposed by the Ref. [55] (Maciula and Szczurek), and also our calculation in the k t -factorization using the MRW UPDFs model of its UPDFs. Finally, it should be clarified that although the collinear framework can describe experimental data well with respect to the k t and (z, k t )-factorizations, it should be noted that these two frameworks can also exhibit remarkable results of exclusive processes with suitable UPDFs and DUPDFs. For instance, our recent works on three-photon productions [23] and Drell-Yan [24] processes have demonstrated that NLO-MRW and PB UPDFs can yield remarkable results. However, in this work, we aim to provide a one-toone comparison between the k t and (z, k t )-factorizations to better understand the difference between these two frameworks. As we have demonstrated in this work, the k t and (z, k t ) frameworks exhibit relatively similar behavior in all regions except for large rapidity limits. One reason we do not compare our predictions with other UPDFs models is that obtaining DUPDFs with those models is a challenging task and there is no straightforward approach to obtain them. Nevertheless, it is an important goal to generate these DUPDFs and better describe the results. and