High-Energy Factorization for Drell-Yan process in $pp$ and $p{\bar p}$ collisions with new Unintegrated PDFs

The formalism for uniform description of Drell-Yan transverse-momentum spectrum is presented in a framework of High-Energy Factorization, which smoothly interpolates between Collins-Soper-Sterman formalism at $|{\bf q}_T|\ll Q$ and usual Collinear Parton Model at $|{\bf q}_T|\sim Q\ll \sqrt{S}$. The new formula for deriving Unintegrated Parton Distribution Functions(UPDFs) from collinear ones is introduced, which leads to excellent description of the shape of $Z$-boson $|{\bf q}_T|$-spectrum at high energies up to $|{\bf q}_T|/\sqrt{S}\simeq 0.02$. Description of normalized $|{\bf q}_T|$-distributions at low energies is achieved via the fit of non-perturbative parameters of quark UPDFs. Reasonable description of angular distributions of leptons in the dilepton center-of-mass frame is also obtained with new UPDFs.


I. INTRODUCTION
The transverse-momentum (q T ) distribution of Drell-Yan(DY) lepton pairs with large invariant mass Q Λ QCD , produced in hadronic collisions, continues to attract a lot of attention from theorists and experimentalists alike. High-precision data on the |q T |-spectrum of lepton pairs with Q close to the Z-boson mass had been obtained very recently by ATLAS Collaboration [1] in pp-collisions with highest energy achieved so far, √ S = 13 TeV. Complimentary set of data on transverse-momentum distribution at lower values of Q had been recently published by PHENIX Collaboration at RHIC Collider with √ S = 200 GeV [2], which partially fills the gap between Drell-Yan data obtained in fixed-target experiments in 1980s and early 1990s [3][4][5] and data obtained at Tevatron [6] and LHC energies.
From the theory side, the description of Drell-Yan |q T |-spectrum at |q T | Q have recently reached maturity, with the achievement [7,8] of Next 3 -to-Leading Logarithmic (N 3 LL) accuracy of the resummation of higher-order perturbative QCD corrections, enhanced by large ln Q 2 /q 2 T , consistently interfaced with non-perturbative effects, important at |q T | ∼ Λ QCD , in the context of the Transverse-Momentum Dependent(TMD)-factorization formalism [9].
At the same time it has been observed [10], that Next-to-Leading Order calculation of the |q T |-spectrum in the Collinear Parton Model(CPM) of QCD can not describe normalization and shape of low-energy Drell-Yan data in the region |q T | > ∼ Q, where ln Q 2 /q 2 T -enhancement of higher-order corrections is not present, and fixed-order predictions should be applicable.
The similar difficulty with the description of transverse-momentum spectrum of identified hadron in Semi-Inclusive Deep-Inelastic Scattering has been found in Ref. [11]. The resummation of threshold effects up to Next-to-Leading Logarithmic Approximation improves the agreement with experimental data only marginally, as it has been shown in Ref. [10]. In our opinion, this phenomenological puzzle is a manifestation of deeper theoretical issue with current formulation of TMD-factorization, which does not provide a unique prescription for the matching between TMD (the so-called W -term) and Collinear-factorization (the Y -term) parts of the calculation at |q T | Q (see e.g. Ref. [12] for detailed discussion) and even lacks QED gauge-invariant definition for the W -term at |q T | ∼ Q, see Refs. [13][14][15].
In the present paper, we approach the problem of uniform description of the |q T |-spectrum of Drell-Yan lepton pairs from a point of view of High-Energy Factorization (HEF), which initially has been introduced as a resummation tool for lnŝ/(−t)-enhanced corrections to the hard-scattering coefficients in Collinear Parton Model [16,17], where invariantsŝ andt refer to the partonic subprocess. Our Parton Reggeization Approach (PRA) is is a version of HEF, based on the Modified Multi-Regge Kinematics (MMRK) approximation for QCD scattering amplitudes. This approximation is accurate both in the Collinear limit, which drives the TMD-factorization and in the High-Energy (Multi-Regge) limitŝ (−t) ∼ q 2 T ∼ Q 2 which is important for Balitsky-Fadin-Kuraev-Lipatov(BFKL) [18][19][20] resummation of lnŝ/(−t)-enhanced effects. This approximation allows us to derive the factorization formula for Drell-Yan cross-section, which is equivalent to the perturbative Collins-Soper-Sterman (CSS) formalism [21] for |q T | Q and accuracy of which at |q T | ∼ Q is expected to increase power-like with decreasing values of |q T |/ √ S and Q/ √ S. Thus, with increasing collision energy we should achieve a uniform description of |q T |-spectrum which does not require any dedicated matching procedure at |q T | Q.
The present paper has the following structure. In the Sec. II we introduce the MMRK approximation and derive factorization formula of PRA for the DY process; in the Sec. III we derive the Unintegrated Parton Distribution Function (UPDF) of PRA; in the Sec. IV we compare our cross-section formula at |q T | Q with the results of CSS formalism up to NLL; in the Sec. V we derive formulas for differential cross-section and squared LO PRA matrix element used in the numerical calculations; in the Sec. VI we compare our predictions with low-energy DY data and perform the fit of non-perturbative parameters of our UPDF; in the Sec. VII we compare our predictions for normalized |q T |-spectra and coefficients parametrizing angular distributions of leptons in the center-of-mass frane of a lepton pair with High-Energy ATLAS [1] and CDF data [6] and in the Sec. VIII we summarize our conclusions.

II. PRA AS HIGH-ENERGY FACTORIZATION AT LEADING POWER
The DY lepton-pair production at leading order in QED coupling constant α proceeds via the exchange of a virtual photon or Z-boson with four-momentum (q), thus the cross-section of this process, differential over invariant mass (Q, q 2 = Q 2 ), rapidity (y) and transversemomentum of the lepton pair, admits a well-known factorization into a convolution of leptonic (L µν ) and hadronic (W µν ) tensors. The latter can be written in the framework of CPM as follows: where is the momentum-density PDF, indices i, j = q,q, g run over parton species, w (ij,aā,CPM) µν is the partonic tensor and index a = V, A(ā = V, A) denote respectively the vector or axial-vector coupling of a vector boson to the quark line in the partonic amplitude (complex-conjugate amplitude). Here and below we use the Sudakov basis-vectors n µ − = 2P µ 1 / √ S and n µ To isolate the x ± -dependence of the cross-section, one introduces the Mellin transform: and then the differential cross-section of a DY lepton pair production via virtual-photon exchange can be written as: where H (ij) and we have introduced dimensionless parameter p 2 = q 2 T /Q 2 T . Note that 0 ≤ p 2 ≤ 1 and p → 0 corresponds to the collinear regime with q 2 T Q 2 , while p = 1 corresponds to the production of the on-shell photon.
Another kinematic limit, in which QCD amplitudes admit simple factorization, is the Regge limit z − → 0, while relation between Q 2 and q 2 T can be arbitrary. In this Multi-Regge Kinematics (MRK), final-state particles are highly-separated in rapidity. Asymptotic expression for QCD amplitudes with quark-exchange int-channel in this limit can be obtained using the formalism of Gauge-Invariant EFT for Multi-Regge processes in QCD [28,29]. For both considered squared amplitudes, this asymptotics depicted diagrammatically in the left panel of the Fig. 1, can be written as: wherek = k µ γ µ , q 2 = q − k 1 ,q 2 = k 2 − k 3 , P µν (q) = −g µν + q µ q ν /Q 2 is the polarization sum for off-shell photon, factor 1/2 corresponds to the averaging over helicities of initial-state quark, Fadin-Sherman scattering vertices [30,31] are: and factorsŜ qi (k 2 ,q 2 ) correspond to the lower part of diagrams in the Fig. 1, with the following general expressions forŜ (±) qi : where factors 1/2 correspond to the averaging over helicities of initial-state quark or gluon, Dirac projectorsP ± =n ∓n± /4 are required by EFT Feynman rules [29] and P µν (k) = −g µν + (k µ n ν + k ν n µ )/(kn) − k µ k ν n 2 /(kn) 2 is the gluon polarization sum in general axial gauge. Note, that due to the structure of vertices (7) and conditions k 2 2 = 0 or (q 2 −k 2 ) 2 = 0, the splitting-factorsŜ qi are invariant w.r.t. the choice of gauge-vector n µ . In the Regge limit z − 1, light-cone components ofq 2 obey the hierarchy:q + 2 q − 2 = z − k − 2 and "small"q + 2 -component is usually neglected in the simplification ofŜ qi (k 2 ,q 2 ). However, this kinematic approximation is not necessary in the case of amplitudes with quark exchange int-channel, because relaxing it does not violate gauge-invariance of the splitting-factors. One can recover theq + 2 momentum component form the on-shell condition (k 2 −q 2 ) 2 = 0:q where we take into account that q 2 T 2 = k 2 T 3 = q 2 T and one finds thatq 2 2 = −q 2 T 2 /(1 − z − ). Substituting the latter approximation for q 2 into Eqns. (8) and (9) one obtains: for both casesŜ (−) qq andŜ (−) qg . Substituting this result into Eq. (6) and calculating the trace: one obtains the Modified MRK Approximation (MMRK) for the considered squared amplitudes: . Note that in the MMRK approximation, we have neglected theq + 2 light-cone component in the "hard process" (virtual photon production vertex in the left panel of Fig. 1), while keeping it in the calculation of theŜ (−) qi . This approximation is more general than Eq. (5), because it is accurate in two limits: q 2 T Q 2 for any z − , and z − 1 for any hierarchy between Q 2 and q 2 T . The MMRK analog of Eq. (4) reads: [32][33][34][35]. The θ-function in the last equation defines the region of applicability of MMRKapproximation to be only the case when rapidity of a virtual photon is larger than the rapidity of a quark or gluon. Indeed, the rapidity of a photon is y γ * = ln(q + /q − )/2 = ln(Q T /(k − 2 z − )), while the rapidity of a final-state parton is y 3 = ln(k + 3 /k − 3 )/2 = ln(|q T |/(k − 2 (1 − z − ))), hence the condition z − < ∆(q 2 T , Q 2 T ) is equivalent to y γ * > y 3 . In the opposite case, the "û-channel" MMRK approximation should be used, which is obtained from approximation above by the replacement z + ↔ z − .
The idea behind MMRK-approximation is not new. It was first successfully applied in in the collinear (dashed curves) and MMRK (solid curves) approximations to the corresponding exact results, obtained with the use of Eqns. (3) and (4), as functions of p = the High-Energy Jets approach [36,37] where a goodt-channel-factorized approximation for QCD amplitudes with emissions of multiple additional partons has been constructed via relaxing of some kinematic constraints in corresponding MRK-asymptotic amplitudes, while preserving their QCD gauge-invariance. Later, the TMD-generalizations of usual DGLAP splitting functions describing the splitting of off-shellt-channel partons have been constructed in Refs. [38][39][40] using the same guiding principles. And recently it has been shown in Ref. [41], that the problem of large NLO corrections for gluon-induced observables in HEF can be solved, if the MMRK approximation for QCD amplitudes is used to construct the UPDF evolution equation and corresponding double-counting subtraction terms at NLO.
In the Fig. 2 we compare the functions H (qq) and H (qg) (5) and MMRK (11) approximations with corresponding exact result obtained by substitution of Eq. (4) into the Eq. (3). One can see, that MMRK approximation provides a reasonable estimate for H (ij) larger values of |q T | the error of MMRK-approximation reaches several tens of percent while staying flat all the way up to p = 1. In contrast to this, the error of collinear approximation rapidly increases when p → 1. Thus, using MMRK-approximation, one can construct the expression for the Drell-Yan |q T |-spectrum with effects of initial-state radiation factorized, which will capture the leading-power in x ± -dependence of the cross-section at least up to Thanks tot-channel-factorized nature of MMRK-approximation and processindependence of splitting-factors (8) and (9) one can derive the factorizaiton formula for the contribution of 2 → 3 partonic process with i, j, i , j = q,q, g, to the hadronic tensor in CPM (1). The MMRK approximation for one of such contributions is depicted diagrammatically on the right panel of the Fig. 2 and for general subprocess of the type (12) the MMRK partonic tensor in Eq. (1), integrated over phase-space of momenta k 3 and k 4 with k 2 3,4 = k + 3,4 k − 3,4 − k 2 T 3,4 = 0 can be written as: where we have introduced integrations over light-cone components oft-channel momenta q 1 , is the usual flux-factor of initialstate partons in CPM factorization formula (1). The MMRK approximation for squared amplitude in Eq. (13) reads: [30,31] which we have applied to the case of Drell-Yan process for the first time in Ref. [42]: as follows: with C (a) q = δ aV e q for the photon-quark coupling and Expression (16) for partonic tensor in HEF is free from any gauge ambiguities at q 2 T ∼ Q 2 , since it exactly satisfies the Ward identity q µ w (kl,aā) µν = q ν w (kl,aā) µν = 0, see also [13,14].
To complete the derivation of the HEF factorization formula one substitutes Eq. (14) to Eq. (13), integrates-out k + 3 , k − 4 ,q − 1 ,q + 2 and k T 3,4 using corresponding delta-functions, then introduces momentum-fractions x 1 = q + 1 /P + 1 and x 2 = q − 2 /P − 2 instead of q + 1 and q − 2 and finally substitutes the result for w (ij,aā) µν into Eq. (1) to obtain: where the tree-level unintegrated PDFs (UPDFs) are : The θ-functions in Eq. (19) enforce the rapidity-ordering between particles in the finalstate y 3 > y γ * > y 4 , for our MMRK approximation for the squared amplitude and kinematics to be applicable. The natural choice of rapidity scale µ Y for the case of Drell-Yan process is µ Y ∼ Q T . As it follows from the discussion above, Eq. (19) is accurate in the region For t µ Y the tree-level UPDF contains a collinear divergence ∼ 1/t signaling the break-down of fixed-order perturbation theory for this object.
An important feature of Eq. (18), which is critical in the region q 2 T ∼ Q 2 S is, that flux-factor 2Sx 1 x 2 is used for off-shell initial-state partons with q 2 1 = −q 2 T 1 < 0 and q 2 2 = −q 2 T 2 < 0. This prescription follows from the derivation of Eqns. (18) and (19), presented above. The similar derivation for gluon-induced processes has been given in Sec.
II of our Ref. [44]. We stress again, that this prescription is necessary for consistency of the cross-section formula of High-Energy Factorization with exact QCD results in Regge limits with x + z + 1 and/or x − z − 1, which give a major contribution to the crosssection in the regime q 2 T ∼ Q 2 S. Thus the prescription of Eq. (18) for the flux-factor should be used consistently with the gauge-invariant amplitudes based on the vertex (15) and both of this factors are important for the |q T |-distribution at |q T | ∼ Q. In connection with this we would like to emphasize that only the "off-shell cross-section" formula (56) in recent Ref. [45], is self-consistent at |q T | ∼ Q, while the "on-shell" cross-section formula (47) is applicable only for |q T | Q.

III. UNINTEGRATED PDF WITH EXACT NORMALIZATION
To resolve a divergence problem of Eq. (19) we follow the standard definition of the UPDF in BFKL formalism (see e.g. Eq. (2.4) in the Ref. [16] or Sec. 1 in [46]) and require that: which is equivelent to: with some function T i (t, µ 2 , x) which is usually referred to as Sudakov formfactor, satisfying the boundary conditions T i (t = 0, µ 2 , x) = 0 and T i (t = µ 2 , µ 2 , x) = 1. We will obtain the latter by multiplying Eq. (19) on the formfactor: and asking for exact equivalence of two definitions (21) and (22). Note, that Eq. (22) coincides with Eq.
Taking the derivative in Eq. (21) with the help of the following from of DGLAP equations for PDFs: which in the limit δ 0 → 0 is exactly equivalent to the usual form of DGLAP equations with (+)-distributions, one obtains: To make contact with Eq. (22), one inserts the identity: into the z-integrands in Eq. (24). Then each integral over z can be split in two terms with integrations over regions x ≤ z ≤ ∆(t, µ 2 ) and ∆(t, µ 2 ) < z ≤ 1 − δ 0 (assuming that ∆(t, µ 2 ) < 1 − δ 0 ) and after reshuffling of some terms, one obtains: In the first line of this equation we have got exactly Eq.(22), therefore we have to put to zero the expression in curly brackets in Eq. (25), which leads to a differential equation for T i (t, µ 2 , x). Another important observation is, that one can safely put δ 0 = 0 in Eq. (25).
Indeed, if i = j then there is no singularity in P ij (z) at z = 1 so integral just converges, while if i = j, then singularity at z = 1 cancels between two terms in the inner-most circular brackets, so integral over z is convergent for δ 0 = 0 anyway.
The solution of obtained differential equation for Sudakov formfactor, satisfying boundary condition T i (t = µ 2 , µ 2 , x) = 1 has the form: with We have written these formulas in the Ref. [47] for the first time, without a detailed derivation. The Sudakov formfactor without the ∆τ i -term in the exponent is similar to the Sudakov formfactor of LO KMRW UPDF of Ref. [35] but with a numerically-important difference that in our MMRK approach, the rapidity-ordering condition is imposed both on quarks and gluons, while in KMRW-approach it is imposed only on gluons. The term proportional to the ratio of PDFs in Eq. (28) is familiar from the expression for parton nonemission porbability in "unitary" Parton Showers [48]. Strictly-speaking, this term makes transformation from PDF to UPDF non-linear w.r.t. the former.
Important property of Eqns. (22), (27), (28) is that they guarantee exact equivalence of definitions (21) and (22) at any order in α s and scheme-choice for DGLAP splitting functions P ij (z) as soon as the PDFsf i (x, µ 2 ) satisfy usual DGLAP equations with the same splitting functions. For alternative ways to ensure the exact normalization condition (20) for KMRW-type UPDF see Ref. [49].

IV. COMPARISON WITH COLLINS-SOPER-STERMAN FROMALISM
For the hadroproduction of Drell-Yan lepton pairs with q 2 T Q 2 the perturbative resummation of higher-order corrections enhanced by ln(Q 2 /q 2 T ) is performed by Collins-Soper-Sterman formula [21]: where x T is a transverse coordinate, conjugated to transverse-momentum q T , C ij are the collinear matching-functions, which are usually taken order-by-order in α s and the resummation is performed by Sudakov formfactor in the x T -space: where functions A and B, corresponding respectively to the resummation of doubly (∝ ln 2 (x 2 T Q 2 )) and single-logarithmic (∝ ln(x 2 T Q 2 )) corrections admit the following perturbative expansion (Eqns. (3.18) and (3.20) in [21]): where we have explicitly shown terms up to Next-to-Leading Logarithmic (NLL) Approximation. The NLL coefficient B in Eq. (32) depends on the resummation scheme, which is defined by parameters C 1,2 in the Ref. [21].
On our momentum-space language, the formfactor (30) The correction ∆τ q is a quantity O(1−∆), so it contributes only beyond NLL-approximation.
Substituting the last result for τ q into Eq.(26) and taking into account, that for t µ 2 : 1 − ∆(t, µ 2 ) t/µ 2 , one obtains in this limit: where we have taken into account, that running-coupling effects in Eq.(26) also contribute only beyond NLL as well as effects of scale-dependence of the PDF in Eq. Eq. (21). So one should consider only the Fourier-transform of a derivative dT qq (q 2 T , µ 2 )/dq 2 T . Taking the Fourier-transform of a t-derivative of Eq. (33) order-by-order in α s , with the help of the relation: 1 where by ellipsis we denote non-logarithmic terms, one obtains: The last result indeed coincides with the square-root of Eq. (30) with coefficients (31) and (32) taken up to NLL-approximation in a scheme with ln C 1 /(2C 2 ) = −γ E . So we conclude, that our resummation scheme is consistent with perturbative part of CSS formalism up to NLL-approximation in the region q 2 T Q 2 where both formalisms apply, thanks to a particular small-t asymptotics of the KMRW cutoff function: 1 − ∆(t, µ 2 ) t/µ 2 .

V. DRELL-YAN LEPTON PAIR PRODUCTION IN PRA
The LO in α s cross-section of p(P 1 ) + p(P 2 ) → l + (p 1 ) + l − (p 2 ) + X-process in PRA is given by: where we have introduced an integration over intermediate momentum q = p 1 + p 2 , parton momenta are given by q µ 1 = x 1 P µ 1 +q µ T 1 and q µ 2 = x 2 P µ 2 +q µ T 2 and PRA squared matrix element |M ij | 2 is a function of scalar products of four-momenta of partons(q 1,2 ), leptons(p 1,2 ) and vectors n + or n − .
In the last line of Eq. (34) one can integrate-out p 2 using the delta-function and then pass to the center-of-mass frame of the lepton pair, to express this integral in terms of spherical angles θ l + and φ l + , parametrizing the direction of lepton momentum in this frame.
In first two lines of Eq.(34) one integrates-out momentum q and momentum-fraction x 2 , using the relation: and replaces d 2 q T 2 → πdq 2 T to finally obtain: where y is the rapidity of the lepton pair, so that x 1,2 = Q T e ±y / √ S and φ 1 is the azimuthal angle of q T 1 , while q T 2 = q T − q T 1 . If four-momenta q µ = (Q T chy, |q T |, 0, Q T shy) µ and q µ 1,2 = ( √ Sx 1,2 /2, q T 1,2 , ± √ Sx 1,2 /2) µ are given in the pp center-of-mass frame, then fourmomenta of leptons can be expressed using covariant relations: with the help of following expressions for pp center-of-mass frame components of unit vectors X µ , Y µ and Z µ of the Collins-Soper frame [50]: The squared PRA amplitude of the LO in α s partonic subprocess: where by Q(Q) we denote Reggeized quark(anti-quark) is given by: where the first, second and third terms in curly brackets correcpond to the squared photon, Z-boson exchange diagrams and Z * γ * -interference respectively. In our numerical calculations we have taken m l = 0, however here we write-down formulas for m l = 0 for generality.
After taking all traces and index-contractions, the squared amplitude can be simplified as follows: where C (q) Zl − is the product of quark and lepton coupling-factors (17), while If instead of the process (36) one considers the processQ(q 1 ) + Q(q 2 ) → l + + l − , the overall sign of the Lorentz-structure M 2 V V,AA should be flipped.

DATA
Perturbative expression (22) does not define the UPDF for all values of t, since for t < Λ 2 QCD , the integral in Eq. (26) is ill-defined due to Landau pole of α s (t ). Similar problem arises also in TMD-factorization [9] in a form of non-perturbative ambiguity of the rapidity-evolution kernel [7]. Analogously with the Ref. [35] we define the UPDF for t < t 0 = 1 GeV as: where parameters A, t 1 and α are determined from the requirements of normalisation of UPDF (20), it's continuity, smoothness in the point t = t 0 and positivity of α as follows: The UPDF defined for all values of t as described above we call the shower-part of the UPDF. Physically it is determined by perturbative dynamics of QCD for t > t 0 and non-perturbative properties of QCD vacuum for t < t 0 [52,53]. To take into account nonperturbative intrinsic motion of partons inside hadron we convolute the shower part on UPDF with phenomenological intrinsic transverse-momentum distribution, which we take in the x-independent Gaussian from: The non-perturbative parameters σ T i will be determined below from a fit of experimental data with fine binning in |q T | < 1 GeV region, but from physical interpretation given above one expects to find σ T i ∼ Λ QCD 1 GeV.
We perform the global fit of parameters σ T i , using experimental data on Drell-Yan lepton pair production at √ S ≤ 200 GeV, summarized in the Tab. I. To obtain the shower part of the UPDF with the procedure described above, we use the LO PDF set MSTW-2008 [54] and formulas (22) Since we do not expect our formalism to describe overall normalization of the data, due to the lack of complete NLO loop corrections, we perform the comparison with normalized distributions (1/σ)dσ/dq 2 T . To this end we multiply each experimental spectrum by constant, q T −independent factor, obtained via the summation of central-values of cross-section in all data-bins. In the Tab. I we show the obtained ratios of experimental total cross-sections and our theoretical results. We also present the uncertainties on this K-factors due to the scalevariation in theoretical predictions (with experimental cross-sections fixed at their central values) and due to experimental uncertainties (divided by central theoretical predictions).
Although in principle, parameters σ T i could be flavor-dependent, due to a large theoretical uncertainty of our LO analysis we have not found any significant improvement in the fit quality from taking different values of σ T i for different falvors or for sea vs. valence quarks.
Therefore we present only the results with all σ T i taken to be equal to the same constant σ T , for which we have found σ binning at |q T | < 1 GeV, the non-perturbative part of the UPDF have negligible effect on our predictions presented in this section. Of course, the usage of NLO UPDF without complete NLO corrections to the PRA hard-scattering coefficient (38) is not fully consistent, however we expect that at least for |q T | < ∼ Q the effect of NLO corrections in PRA on the shape of the distribution will be negligible and NLO correction will affect mostly the overall normalization of the spectrum.
In the Fig. 6 we compare our predictions for normalized distribution (1/σ)dσ/d|q T | with At higher transverse momenta the power-corrections w.r.t. x ± become important, which can be taken into account only by the complete NLO calculation in PRA. In conjunction with this results we can also point towards the recent study [59], where UPDFs defined by Eq. (21) with x-independent Sudakov formfactor and gauge-invariant Matrix Elements with off-shell initial-state partons derived in a formalism, which is equivalent to ours [60,61], had been successfully used to describe the Z-boson production in proton-lead collisions.
In the Fig. 8 we compare our quark UPDFs with UPDFs obtained in the receltly-proposed Parton-Branching(PB) method [62,63]. The latter UPDFs can be obtained from the TMDlib package [64]. The UPDFs in PB-method are derived as Monte-Carlo solution of a system of evolution equations constructed in such a way, that q T -integrated UPDF satisfies usual DGLAP equations, while transverse-momentum dependence of UPDF is essentially determined from the ambiguity in definition of "non-resolved" parton branchings by means of a suitable cutoff function and several scale-choices in the definition of Sudakov formfactor and branching probability. From the Fig. 8   estimate the cross-section. To overcome this problem, authors of recent Ref. [67] attempt to match the LO calculation with PB UPDFs with NLO QCD corrections obtained via the standard implementation of the MC@NLO method [68]. In such a way, satisfactory description of shapes and normalization of low-energy Drell-Yan data, as well as ATLAS data on Z-boson q T -spectrum for |q T | < 10 GeV has been obtained in Ref. [67]. However the standard MC@NLO method is not designed to properly take into account the off-shell initial state partons and hence it's application together with UPDFs is hard to justify. The consistent formalism of NLO calculations with off-shell initial-state partons is currently under development [41,47,[69][70][71][72]. Moreover, taking into account the power-like tail of UPDF, together with effects of quark Reggeization, allows one to extend the range of applicability of High-Energy Factorization for Z-boson production at √ S = 13 TeV all the way up to |q T | < ∼ 200 GeV as we have shown above.
In the recent Ref. [45] the gluon and quark UPDFs have been constructed as a solutions of CCFM-Kwieciński evolution equation, and consistency of this formalism with CSS-formalism up to NLL-approximation has been demonstrated. A good description of shapes of Drell-Yan |q T |-spectra at low energies, similar to ours (Figs. [3][4][5], has been obtained in Ref. [45]. It is interesting to note, that evolution equations in this formalism are very similar to the PB evolution equations, while the logic to obtain them is different. It is tempting to suggest, that closed-form solution of CCFM-K or PB evolution equations could be found in terms of underlying collinear PDFs, analogous to our Eqns. (22), (26), (27) and (28).
Finally, we shall discuss how our approach describes the angular distribution of leptons in the rest-frame of the lepton-pair, which can be parametrised as follows: where angular coefficients A i are functions of Q 2 , q 2 T and y, realted with the polarization density-matrix of the intermidiate vector-boson. Thus the study of transverse-momentum dependence of angular coefficients will allow us to check whether the spin structure of our MMRK amplitudes can reasonably approximate the spin structure of exact QCD amplitudes.
To obtain theoretical predictions for angular coefficients with the help of our masterformula for the cross-section (35), we use the same harmonic-projectors method which has been used to obtain theoretical predictions in the Ref. [51], see Eqns. (4) and (5) of the the Sec. 2 in this reference.
In the Fig. 9 we compare our theoretical predictions for angular coefficients with experimental data obtained by ATLAS Collaboration [51] in pp-collisions with √ S = 8 TeV in the range of lepton-pair invariant masses 80 < Q < 100 GeV and lepton-pair rapidities |y| < 2.
The results on coefficients A 0 and A 2 are especially interesting since in the CPM the Lam-Tung relation [73,74]: A 0 = A 2 holds up to NLO in α s and is violated only by NNLO effects. The coefficients A 1 , A 3 and A 4 are nonzero only if parity-violating couplings of Z-boson are present, see the right panel of Fig. 9, while coefficients A 5 , A 6 and A 7 further require effects beyond NLO of CPM to be taken into account and these coefficients are zero in the LO of PRA. From the right panel of Fig. 9 one observes, that q T -dependence of A 3 is nicely described in the LO of PRA, while A 1 and A 4 come-out to be of the right order of magnitude and roughly of a correct shape. Nevertheless, it is clear, that LO of PRA is not capable to correctly capture this subtle details of a spin structure of production amplitude, and full NLO corrections in PRA, which will exactly take into account emission of an additional parton from the hard process, are necessary to quantitatively predict all angular coefficients.

VIII. CONCLUSIONS AND OUTLOOK
In the present paper we have introduced a new prescription to obtain UPDF from LO and NLO collinear PDFs and to define it's non-perturbative ambiguity. This UPDF, together with QCD and QED gauge-invariant matrix elements already in the LO in α s provide an excellent description of shapes of transverse-momentum distributions of DY lepton-pairs in pp and pp-collisions at low and high collision energies, in the region Q T / √ S 1. Qualitative description of transverse-momentum dependence of angular coefficients of lepton distribution in the rest frame of the lepton pair is also achieved for Q M Z . To describe the normalization of |q T |-distribution and extend the formalism to higher values of transverse momenta it is necessary to go beyond LO in α s for the coefficient function in our approach.