Relation between the parton branching approach and Catani-Ciafaloni-Fiorani-Marchesini evolution

We consider the associated production of electroweak gauge bosons and charm or beauty quark jets at the LHC using the kT -factorization framework. We apply the transverse momentum dependent (TMD) parton distributions in a proton obtained from the parton branching (PB) method as well as from the CataniCiafaloni-Fiorani-Marchesini (CCFM) evolution equation. For the PB approach, our prescription merges the standard leading order OðααsÞ kT-factorization calculations with several tree-level next-to-leading order OðααsÞ off-shell production amplitudes. For the CCFM scenario, our consideration is based on the OðααsÞ off-shell gluon-gluon fusion subprocess g g → Z=W þQQ̄ and some OðααsÞ subprocesses involving quark interactions, taken into account in conventional (collinear) QCD factorization. We find that the W þ c and Z þ b cross sections, calculated within the PB and CCFM-based schemes with the proper choice of leading and next-to-leading subprocesses, are in good agreement with each other, thus establishing a correspondence between these two scenarios. A comparison with the latest LHC experimental data is given and the necessity for the proper off-shell treatment of the production amplitudes in determination of the parameters of the TMD parton density is demonstrated.


I. MOTIVATION
A theoretical description of a number of processes at high energies and large momentum transfer containing multiple hard scales requires so-called transverse momentum dependent (TMD, or unintegrated) parton (quark and gluon) density functions [1]. These quantities encode the nonperturbative information on proton structure, including transverse momentum and polarization degrees of freedom and are related to the physical cross sections via different TMD factorization scenarios. The latter provide the necessary framework to separate hard partonic physics, described with a perturbative QCD expansion, from soft hadronic physics.
In the limit of a fixed hard scale and high energy the k T -factorization [2] (or high-energy factorization [3]) approach is expected to be valid. This approach is mainly based on the Balitsky-Fadin-Kuraev-Lipatov (BFKL) [4] or Ciafaloni-Catani-Fiorani-Marchesini (CCFM) [5] evolution equations for the TMD gluon densities in a proton.
The BFKL equation resums large logarithmic terms proportional to α n s ln n s ∼ α n s ln n 1=x, important at high energies s (or, equivalently, at small x). The CCFM equation takes into account additional terms proportional to α n s ln n 1= ð1−xÞ and is valid at both low and large x. There are also scenarios to evaluate the TMD parton densities based on the conventional Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) [6] evolution equations, namely the Kimber-Martin-Ryskin (KMR) prescription [7] and recently proposed parton branching (PB) approach [8,9]. The KMR approach, currently explored [10] at next-to-leading order (NLO), is a formalism to construct the TMD parton densities from well-known conventional (collinear) ones under the key assumption that the transverse momentum dependence of the parton distributions enters only at the last evolution step. The PB approach provides an iterative solution of the DGLAP evolution equations for collinear and TMD parton density functions by making use of the concept of resolvable and nonresolvable branchings and by applying Sudakov form factors to describe the parton evolution from one scale to another without resolvable branching. The splitting kinematics at each branching vertex is described by the DGLAP equations and angular ordering conditions for parton emissions can be applied here instead of the usual DGLAP ordering in virtuality. One of the main advantages of the PB approach is that the free parameters can be determined from fits to inclusive measurements and when later these PB TMD densities are used in a parton shower approach, no further adjustment of parameters is needed (and allowed). 1 A number of phenomenological applications of the CCFM evolution equation is known in the literature (see, e.g., [11][12][13][14][15][16][17] and references therein). Several applications of the PB approach were discussed [18,19] and comparison between the PB and KMR predictions has been recently made [20,21]. However, the correspondence between the CCFM and PB scenarios has not been investigated yet. One of the main goals of this paper is to compare predictions for some QCD processes with the CCFM and PB parton distributions, to find a correspondence between these approaches and to define conditions, at which such correspondence takes place. As the reference processes for the study, we consider here the associated production of electroweak gauge bosons (W and Z) and charm or beauty quark jets at the LHC conditions. These are the so-called "rare" processes which could have never been systematically studied at previous accelerators. We already successfully applied [16] the CCFM-evolved gluon densities to describe the first LHC data [22,23] of the associated Z þ b production at ffiffi ffi s p ¼ 7 TeV. Those calculations were based on the Oðαα 2 s Þ off-shell gluon-gluon fusion subprocess g Ã g Ã → Z þ QQ (where the Z boson further decays into a lepton pair) and several subleading Oðαα 2 s Þ and Oðαα 3 s Þ subprocesses involving quarkantiquark and quark-gluon interactions, taken into account within the conventional (collinear) QCD factorization. Such a scheme allows us to describe LHC experimental data in the whole transverse momentum range (see [16] and Sec. II C below for a more detailed explanation of our calculation scheme). Here we extend the consideration to associated W þ c production, measured [24] by the CMS collaboration for the first time as a function of W decay lepton and/or c-jet rapidities at ffiffi ffi s p ¼ 13 TeV. Thus, in this sense we continue the line of our previous studies [16]. In contrast to the CCFM scenario, in the PB calculations (as being the DGLAP-based ones) one has to include usual leading order (LO) Oðαα s Þ subprocesses properly matched with a number of additional higher-order terms. Below we perform such calculations and the matching procedure following the approach applied recently [20] for c-jet production at the LHC. The comparison between the results obtained within the TMD approaches above could also be a general consistency check for the k T -factorization phenomenology. The outline of the paper is the following. In Sec. II we briefly discuss the CCFM equation and PB approach and describe the basic steps of our calculations. In Sec. III we present the results of our calculations and discussion. Our conclusions are summarized in Sec. IV.

A. CCFM evolution
The CCFM gluon evolution equation resums large logarithms α n s ln n 1=ð1 − xÞ in addition to BFKL ones α n s ln n 1=x and introduces angular ordering of initial emissions to correctly treat gluon coherence effects. In the limit of asymptotic energies, it is almost equivalent to BFKL, but also similar to the DGLAP evolution for large x [5]. In the leading logarithmic approximation, the CCFM equation for TMD gluon density Aðx; k 2 T ; μ 2 Þ with respect to the evolution (factorization) scale μ 2 can be written as where k 0 T ¼ qð1 − zÞ þ k T , scale μ 0 is the soft starting scale of the evolution and z M ¼ 1 − μ 0 =μ 0 is a resolution parameter. 2P gg ðz; k 2 T ; q 2 Þ is the CCFM splitting function: The Sudakov and non-Sudakov form factors read whereᾱ s ¼ 3α s =π, and the z-integral in (4) is finite due to the theta functions [25]. The first term in the CCFM equation, which is the initial TMD gluon density multiplied by the Sudakov form factor, corresponds to the contribution of nonresolvable branchings between the starting scale μ 2 0 and scale μ 2 . The second term describes the details of the QCD evolution expressed by the convolution of the CCFM gluon splitting function with the gluon density and the Sudakov form factor. The theta function introduces the angular ordering condition. The evolution scale μ 2 is defined by the maximum allowed angle for any gluon 1 In standard parton shower event generators, the parton shower parameters have to be determined separately. 2 It was shown in [9] that with this parameter virtual and resolvable branchings are treated consistently as long as z M → 1. emission [5]. A similar equation also can be written [26] for valence quark densities 3 (with replacement of the gluon splitting function by the quark one). Usually, the initial TMD gluon and valence quark distributions are taken as where σ ¼ μ 0 = ffiffi ffi 2 p and q v ðx; μ 2 Þ is the standard (collinear) density function. The parameters of the initial TMD parton distributions can be fitted from the collider data (see, e.g., [28,29]).
The CCFM equation can be solved numerically using the UPDFEVOLV program [30], and the TMD gluon and valence quark densities can be obtained for any x, k 2 T and μ 2 values. The main advantage of this approach is the ease of including higher-order radiative corrections (namely, a part of NLO þ NNLO þ Á Á Á terms corresponding to the initialstate real gluon emissions) even within LO. More details can be found, for example, in review [1].

B. Parton branching approach
The parton branching approach was introduced [8,9] to treat the DGLAP evolution. The method provides an iterative solution of the evolution equations and agrees with the usual methods to solve the DGLAP equations for inclusive distributions at the NLO and NNLO. It allows one to take into account simultaneously soft-gluon emission at z → 1 and transverse momentum q T recoils in the parton branchings along the QCD cascade. The latter leads to a natural determination of the TMD quark and gluon densities. A soft-gluon resolution scale z M is introduced to separate resolvable and nonresolvable emissions, which are treated via the DGLAP splitting functions P ab ðα s ; zÞ and Sudakov form factors, respectively. The PB equations for TMD parton densities read The evolution scale μ 2 can be connected with the angle of emitted parton with respect to the beam direction, which leads to the well-known angular ordering condition, μ ¼ jq T j=ð1 − zÞ. It was shown in [9] that stable results for the collinear distributions are obtained if z M is close to one, which is the case for the PB distributions used here. The initial TMD parton distributions are taken in a factorized form as a product of collinear quark and gluon densities and intrinsic transverse momentum distributions (treated as Gaussian ones [8,9]), where all the parameters can be fitted from the collider data. Unlike the CCFM parton distributions, the PB densities have the strong normalization property: The PB evolution equations can be solved by an iterative Monte-Carlo method, which results in a steep drop of the parton densities at k 2 T > μ 2 . It contrasts the CCFM evolution, where the transverse momentum is allowed to be larger than the scale μ 2 , corresponding to an effective taking into account higher-order contributions. 4 C. Associated W AE =Z + Q production with TMD factorization To calculate total and differential cross sections of associated gauge bosons and heavy quark jet production within the CCFM-based approach, we strictly follow the scheme [16]. In this scheme, the leading contribution (in the small-x region, where the gluon density dominates over the quark distributions) comes from the Oðαα 2 s Þ off-shell gluon-gluon fusion subprocess In addition to off-shell gluon-gluon fusion, one can take into account several Oðαα 2 s Þ subprocesses involving quarks in the initial state: for Z þ b production and for W − þ c production. Subprocesses for W þ þc production can be obtained via charge conjugation. 5 The quarkinduced diagrams may become important at very large transverse momenta (or, respectively, at large x, which is needed to produce large p T events), where the quarks are less suppressed or can even dominate over the gluon density. Following [16], the contributions from subprocesses (12), (13), (15) and (16) are taken into account using the collinear DGLAP-based factorization scheme, which provides better theoretical grounds in the large x region. 6 So, we consider a combination of two techniques with each of them being used at the kinematic conditions where it is best suitable. As it was already mentioned above, the main advantage of this approach is that a large piece of higher order (collinear) corrections, even beyond the NLO order, is included into the calculations via the off-shell gluongluon fusion subprocess supplemented with the CCFM gluon dynamics. Such calculations, although being formally only at leading logarithmic accuracy at low and moderate transverse momenta, result in reasonable pQCD predictions, as it was argued in [32]. At large transverse momenta our scheme employs usual tree-level NLO pQCD evaluations. We note that the contributions from the offshell Oðαα s Þ subprocesses, namely, in the CCFM scheme are covered by gluon-fusion subprocesses (10) or (11) and therefore not taken into account to avoid the double counting. In contrast, in the DGLAPbased PB approach one has to take into account the offshell Oðαα s Þ subprocesses (18) or (19) and properly match them with the higher-order contributions (12)- (14) or (15)- (17), thus keeping the tree-level NLO accuracy. The details of the matching procedure are discussed below (see Sec. III A). Note that the effects of the initial state parton showers are already included into these PB predictions.
Technically, to calculate the cross sections of processes under consideration we have to convolute the relevant partonic cross sections (related with the off-shell production amplitudes) and TMD parton densities in a proton, according to the k T -factorization prescription: where x 1 and x 2 are the longitudinal momentum fractions of the initial off-shell partons and k 2 1T and k 2 2T are their transverse momenta. In the case of gluon-gluon fusion subprocesses (10) and/or (11) we have × dp 2 1T dp 2 2T dy 1 dy 2 dy dϕ 1 2π where ffiffi ffi s p is the total energy, p 1T , p 2T , y 1 , y 2 , ψ 1 and ψ 2 are the transverse momenta, rapidities and azimuthal angles of outgoing quarks, ϕ 1 and ϕ 2 are the azimuthal angles of initial off-shell gluons and y is the rapidity of produced gauge boson. For the quark-antiquark annihilation subprocess (12) one can write × dp 2 2T dy 1 dy 2 dy dϕ 1 2π Similar expressions for other 2 → 3 quark involved subprocesses (13)-(17) can be easily obtained from the factorization formula (20). For the 2 → 2 subprocess (18) we have × dp 2 1T dy 1 dy dϕ 1 2π A similar formula can be easily obtained for remaining subprocess (19). The gauge-invariant off-shell production amplitudes for gluon-gluon fusion subprocesses (10) and (11) were calculated earlier [33,34] and implemented into the Monte-Carlo event generator CASCADE [32] and newly developed parton-level Monte-Carlo event generator PEGASUS [35]. The off-shell amplitudes for quark induced subprocesses (12)- (14) and (15)- (17) can be derived in the framework of the Reggeized parton approach [36]. One can also use Britto-Cachazo-Feng-Witten recursion for offshell gluons [37] and the method of auxiliary quarks for off-shell quarks [38], implemented in the Monte-Carlo generator KATIE [39]. In this study, to calculate the contributions from (12)- (14) and (15)- (17) subprocesses in the PB scheme we used the KATIE tool. To calculate the contributions from the quark-induced subprocesses (12)- (17) in the CCFM scenario, we apply the on-shell limit of the expressions listed above.
In the present paper we compare the CCFM and PB results obtained with JH'2013 set 1 [28] and PB-NLO-HERAI+II-2018 set 2 [19] TMD parton densities in a proton. 7 The main motivation for our choice is that the input parameters of both these distributions were obtained in exactly the same way: from the best description of precision deep inelastic scattering data on the proton structure functions F 2 with exactly the same angular ordering conditions (see [19,28] for more information). For the conventional quark and gluon densities we use the MMHT'2014 (LO) set [41].

III. NUMERICAL RESULTS
Before we show the results of our calculations, we list the input parameters. We use two-loop running strong coupling formula with n f ¼ 5 massless quark flavors and take Λ QCD ¼ 200 MeV in the CCFM case and Λ QCD ¼ 118 MeV for PB distributions [19,28]. The QED running coupling is applied with αðm 2 Z Þ ¼ 1=128. The electroweak bosons masses were taken as m W ¼ 80.385 GeV and m Z ¼ 91.188 GeV [42]. As it is often done, we keep the factorization and renormalization scales to be equal to the gauge boson mass. However, in the CCFM scheme we use a different value for factorization scale μ 2 F ¼ŝ þ Q 2 T , whereŝ and Q T are the energy of the scattering subprocess and transverse momentum of the incoming off-shell gluon pair. The definition of μ F is unusual and dictated by the CCFM evolution algorithm [28].

A. Matching O(αα s ) and O(αα 2 s ) terms in the PB approach
As it was already mentioned above, we supplemented the LO contributions (18) or (19) with off-shell partons in the PB calculations by the tree-level Oðα 2 S Þ corrections (13) and (14) or (16) and (17) from the emission of one additional parton. However, as it is well known, a problem of possible double counting can occur when mixing different final states. Let us consider the Z þ b production (of course, the same arguments apply for W þ c case). Here, the off-shell subprocess (18) partially includes contributions from (13) and (14) due to initial state parton radiation, which can result in substantial double counting if these contributions are summed up. To avoid this double counting we limit the integration over the transverse momenta of the incoming off-shell quark and gluon in the factorization formula (20) for the LO subprocess (18) from above with some value k cut T , so jk T1 j < k cut T and jk T2 j < k cut T . Thus, one removes jets, originating from the initial state radiation and being harder than the initial partons. The latter, however, could be covered by the subprocesses (13) and (14), if we choose there only the events with final gluons and light quarks, having transverse momenta p T larger, that the cut scale k cut T . In this way, therefore, we can almost avoid the double counting region. 8 Of course, the value of k cut T is not universal but is depending on the process. In order to determine k cut T we have calculated the differential cross sections of the LO subprocesses (18) or (19) as a function of initial gluon transverse momentum and differential cross sections of the Oðα 2 S Þ subprocesses (13) or (14) as a function of the produced gluon transverse momentum p T , see Fig. 1. These calculations were performed in the fiducial kinematical region covered by the ATLAS [23] and CMS [24] experiments (see below). So, fixing the k cut T at some value would mean that we keep the contribution from the Oðα S Þ subprocesses lying to the left from the vertical line with k T ¼ k cut T and complement it with the contribution from the Oðα 2 S Þ subprocesses right to the vertical line. The resulting matched p T distribution will have a steplike discontinuous behavior at k T ¼ k cut T . A reasonable choice for the k cut T would be then the one with the step being small. We find that this can be achieved with k cut T ≃ 30 GeV for associated Z þ b-jet production and k cut T ≃ 15 GeV for W þ c-jet production. As one can see from Fig. 1, this choice will lead to the continuous merged transverse momentum distributions.
To investigate the dependence of PB predictions on the k cut T value in more detail, we calculated the ratios of fiducial cross sections σ PB ðZþbÞ=σ CCFM ðZþbÞ and σ PB ðW þ cÞ= σ CCFM ðW þ cÞ as a function of k cut T (note that the denominators in these ratios do not depend on the k cut T ). Our results are shown in Fig. 2. We find that the matched PB predictions are rather stable with variation in k cut T : both 7 A comprehensive collection of the available TMD parton densities can be found in the TMDLIB package [40]. 8 One can estimate the double counting contribution considering the subprocess (14) with the condition jk T1;2 j < k cut T and reconstructing the parton shower with CASCADE [32]. Then one can choose only the events containing a hard outcoming gluon with the transverse momentum p T > k cut T . That would contain events, which are supposed to be covered with the subprocess (14), thus making the double counting. We have evaluated this cross section with k cut T ¼ 30 GeV and found that it is of order ∼5% of the total contribution from the subprocesses (14) and (18) in all the p Z T and y Z subregions of the fiducial region.
the Z þ b-jet and W þ c-jet cross sections change less than 5% if k cut T ≥ 10 GeV or k cut T ≥ 20 GeV, respectively. The dependence on k cut T is smaller than the scale uncertainties of our calculations (see estimation below), so we employ the matching value k cut T ¼ 30 GeV for Z þ b calculations and k cut T ¼ 15 GeV for W þ c production in our numerical calculations.
One can see that with the appropriate choice of k cut T , as discussed above, the fiducial cross sections calculated in the PB approach are very close to the ones obtained in the CCFM scheme. The correspondence between these approaches is investigated in detail in the next section.

B. Comparison with the LHC data
We are now in a position to present the results of our simulations and to confront them with the latest LHC data.
The measurements of the associated production of Z bosons and beauty jets have been carried out by the ATLAS [22] and CMS [23] collaborations and refer to the following categories: Z bosons produced in association with one beauty jet, Z bosons produced in association with two beauty jets, Z bosons associated with any number of b-jets and Z bosons produced in association with explicitly reconstructed b-hadrons. The data on the associated production of W bosons and one charmed jet were reported by the CMS collaboration very recently for the first time [24].
In the present study we concentrate on the production of gauge bosons associated with one heavy quark jet.
The ATLAS collaboration has collected the data on Z þ b-jet production at the center-of-mass energy ffiffi ffi s p ¼ 7 TeV [22]. Both leptons originating from the Z boson decay are required to have p l T > 20 GeV and jη l j < 2.4, the lepton pair invariant mass lies in the interval 76 < M ll < 106 GeV, the beauty jets are required to have p b T > 20 GeV and jη b j < 2.4. The measurement of W þ cjet production at LHC was made by the CMS collaboration [24] at ffiffi ffi s p ¼ 13 TeV and the fiducial region was defined with the following cuts: the transverse momentum of the cquark p c T > 5 GeV, transverse momentum and pseudorapidity of the muon originating from W decay should be p μ T > 26 GeV and η μ < 2.4. The transverse mass of the W boson should be m T > 50 GeV.
We start from Z þ b production. In Fig. 3 we present the predictions for the Z boson rapidity and transverse momentum distributions in comparison with the ATLAS [22] data. The solid and dashed histograms correspond to the results obtained with the JH'2013 set 1 and PB-NLO-HERAI+II-2018 set 2 parton densities. The shaded band represents the scale uncertainties of our PB-based calculations, which have been estimated by varying the scales μ R and μ F by a factor of 2 around their default values. 9 For comparison, we also show the conventional (collinear) NLO pQCD predictions, taken from [22], calculated with MCFM generator [43]. 10 One can see that the Z boson rapidity distribution shows almost perfect agreement between the CCFM and PB approaches, demonstrating the consistency between the CCFM and PB approaches. The cross sections are lower than the ones obtained in the collinear approach. However, the k T -factorization based calculations are in better agreement with the data, though slightly overestimating the ATLAS data in the central region. This overestimation is, however, covered by the uncertainties of our calculations. More information can be obtained with the p Z T -distributions. The CCFM and PB-based calculations give almost the same very good description of the ATLAS data at small p Z T , while at p Z T ≳ 100 GeV the PB-based cross section lies significantly higher than the CCFM one, and is in better agreement with the data. The reason for that is a more accurate treatment of quark-initiated subprocesses in the PB approach, including also contributions from subprocesses (14) and (17). We would like to note that both approaches describe the ATLAS data generally better than the NLO collinear factorization results, especially at low  [22]. 9 The uncertainties connected with the determination of PB TMD parameters are typically much less than the estimated scale uncertainties and are not taken into account in this work. The uncertainties connected with the different PB sets are out of our present consideration. 10 Note that the MCFM predictions were obtained at a different scale m Z T ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m 2 Z þ p Z2 T p and corrected for the effects of QED final-state radiation, hadronization, underlying events and multiparticle interactions. p Z T . The scale uncertainties, estimated for the PB scheme, are even less than the ones of collinear NLO predictions. Now we turn to the W þ c-jet production. Our results are shown in Fig. 4, where we plot the decay muon rapidity distributions for both W þ þc and W − þ c events measured by the CMS collaboration [24]. As in the Z þ b production, one can see that the results obtained with the CCFM and PB approaches agree very well with each other. This confirms once again the consistency between these approaches. However, our predictions are lower than the ones obtained in the collinear QCD factorization using the MCFM tool, which is in contrast to the Z þ b-jet production. The results of MCFM calculations agree well with the CMS data. To explain the reason of the observed underestimation in the TMD predictions, let us consider the relevant differential cross sections as functions of longitudinal momentum fractions and transverse momenta of incoming partons. We show these cross sections on Figs. 5 and 6. As one can see, the W þ c-jet production is dominated by smaller x and broader k T regions in comparison with the Z þ b-jet production case. The corresponding off-shell production amplitudes are known to be suppressed in the large k T domain [44]. We demonstrate this effect in Fig. 7 (left), where the ratios of the reduced cross sections of (18) and (19) calculated with off-shell and on-shell production amplitudes are presented. We find that at k T ∼ 20 GeV the off-shell amplitude becomes greatly suppressed. However, a large part of the W þ c events comes from the region of k T ∼ 20 GeV, as one can see from Fig. 6. Since the considered TMD parton densities were fitted with on-shell matrix elements (see [19,28]), the suppression results in the drop of the total cross section, thus leading to the observed underestimation. The observed flat behavior at relatively low k T is connected with the kinematical cuts applied in the CMS analysis [24]. To investigate it in more detail we repeat the calculations with different cuts on the transverse momentum of the produced c-quark without any other restrictions on the phase space. One can see that with increasing the p cut T the plateau continues until a larger value of k T ∼ p cut T . In the case of a small value of p cut T we obtain a steep behavior starting from practically zero k T . A similar observation was made in [44] for charm and beauty quark photoproduction at HERA, where the heavy quark mass played the role of p cut T . Thus, we conclude that to describe the overall normalization of W þ c-jet production [24] the appropriate fit of the parameters of considered TMD parton distributions with proper off-shell treatment of production amplitudes is needed.
Finally, we also make a prediction for the p W T -distribution in the W þ c production case (Fig. 8). One can see that the CCFM and PB approaches result in different shapes of the distribution; however, the position of the peak remains the same. Like in the case of Z þ b production, the difference in the shapes can be explained by a more accurate treatment of the quark-initiated subprocesses within the PB calculations.

IV. CONCLUSION
We have studied the associated production of Z and W bosons and charm or beauty quark jets at the LHC conditions using the TMD factorization framework. We have applied the TMD parton distributions in a proton obtained from the recent PB method as well as from the CCFM evolution equation. For the PB approach, our prescription merges the Oðαα s Þ calculations with several tree-level Oðαα 2 s Þ off-shell production amplitudes. For the CCFM scenario, our consideration is based on the Oðαα 2 s Þ off-shell gluon-gluon fusion subprocess g Ã g Ã → Z 0 =W AE QQ and some subleading Oðαα 2 s Þ subprocesses involving quark interactions, taken into account in conventional QCD factorization. We have found that the W þ c and Z þ b cross sections, calculated within the PB and CCFM-based schemes with the proper choice of leading and next-to-leading subprocesses in the k T -factorization, are in good agreement with each other. We have demonstrated the necessity for the proper off-shell treatment of the production amplitudes in determination of the parameters of the TMD parton densities in a proton.