On relation between Parton Branching Approach and CCFM 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 Catani-Ciafaloni-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(αα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 gg → Z/W + QQ̄ and some subleading O(αα2 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. PACS number(s): 12.38.-t, 12.38.Bx, 13.85.Qk, 14.65.-q, 14.70.Hp ∗e-mail: malyshev@theory.sinp.msu.ru


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 non-perturbative 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 highenergy 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 usual DGLAP ordering in virtuality. One of the main advantages of the PB approach is that the TMD parton densities (and all corresponding non-perturbative parameters) can be fitted to experimental data, so that the theoretical predictions, where the parton shower effects are already taken into account, can be obtained with no further free parameters 1 .
A number of phenomenological applications of the CCFM evolution equation is known in the literature (see, for example, [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 been not 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 succesfully applied [16] the CCFM-evolved gluon densities to describe first LHC data [22,23] of the associated Z + b production at √ s = 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 quark-antiquark 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. 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 √ s = 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 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 be also a general consistency check for the k T -factorization phenomenology. The outline of the paper is the following. In Section 2 we briefly discuss the CCFM equation and PB approach and describe the basic steps of our calculations. In Section 3 we present the results of our calculations and discussion. Our conclusions are summarised in Section 4.

Theoretical framework 2.1 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 logarthmic 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 ′ T = q(1 − z) + k T andP gg (z, k 2 T , q 2 ) is the CCFM splitting function: The Sudakov and non-Sudakov form factors read: whereᾱ s = 3α s /π. 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 non-resolvable 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 emission [5]. A similar equation also can be written [26] for valence quark densities 2 (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 / √ 2 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, for example, [26,28]).
The CCFM equation can be solved numerically using the updfevolv program [29], 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 initial-state real gluon emissions) even within LO. More details can be found, for example, in review [1].

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 non-resolvable 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: where The evolution scale µ 2 can be connected with the angle of emitted parton with respect to the beam direction, that leads to the well-known angular ordering condition, µ = |q T |/(1 − z). The dependence on the z M falls out when this angular ordering condition is applied and z M is large enough. 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, that 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 3 .

Associated W ± /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 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 4 . The quark-induced 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 5 . So, we consider a combination of two techniques with each of them being used at the kinematic conditions where it is best suitable. We note that the contributions from the off-shell 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 DGLAP-based PB approach one has to take into account the off-shell O(αα s ) subprocesses (18) or (19) and properly match them with the higher-order contributions (12) - (14) or (15) - (17), respectively. The details of the matching procedure are discussed below (see Sec. 3.1). According to k T -factorization prescription, 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: 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. The gauge-invariant off-shell production amplitudes for gluon-gluon fusion subprocesses (10) and (11) were calculated earlier [31,32] and implemented into the Monte-Carlo event generator cascade [33] and newly developed parton-level Monte-Carlo event generator pegasus [34]. The off-shell amplitudes for quark induced subprocesses (12) -(14) and (15) -(17) can be derived in the framework of the reggeized parton approach [35]. One can also use Britto-Cachazo-Feng-Witten (BCFW) recursion for off-shell gluons [36] and method of auxilliary quarks for off-shell quarks [37], implemented in the Monte-Carlo generator katie [38]. In this study, to calculate the contributions from (12) -(14) and (15) -(17) subprocesses in the PB scheme we used the katie tool.
In the present paper we compare the CCFM and PB results obtained with JH'2013 set 1 [26] and PB-NLO-HERAI+II-2018 set 2 [19] TMD parton densities in a proton 6 . 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 DIS data on the proton structure functions F 2 with exactly the same angular ordering conditions (see [19,26] for more information). For the conventional quark and gluon densities we use the MMHT'2014 (LO) set [40].

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 flavours and take Λ QCD = 200 MeV in CCFM case and Λ QCD = 118 MeV for PB distributions [19,26]. 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 [41]. 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 [26].
3.1 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, that 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 |k T 1 | < k cut T and |k T 2 | < 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.
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 step-like discontinuous behaviour 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 details 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 the Z + b-jet and W + cjet 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.

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 centerof-mass energy √ s = 7 TeV [22]. Both leptons originating from the Z boson decay are required to have p l T > 20 GeV and |η l | < 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 |η b | < 2.4. The measurement of W + c-jet production at LHC was made by the CMS Collaboration [24] at √ s = 13 TeV and the fiducial region was defined with the following cuts: the transverse momentum of the c-quark p c T > 5 GeV, transverse momentum and pseudo-rapidity 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 corresponds 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 7 . For comparison, we also show the conventional (collinear) NLO pQCD predictions, taken from [22], calculated with mcfm generator [42]. One can see that the Z boson rapidity distribution show 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 p Z T . The scale uncertainties, estimated for 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 collinear predictions, which is in contrast to the Z + b-jet production. To explain the reason of the observed underestimation, 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 supressed in the large k T domain [43]. 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,26]), 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 details 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 continuates 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 behaviour starting from practically zero k T . A similar observation was made in [43] 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.

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 ± 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. Thus we have established a correspondence between these two scenarios. We have demonstrated the necessity for the proper off-shell treatment of the production amplitudes in determination of the parameters of the TMD parton dentities in a proton.        densities. M.A.M. was also supported by a grant of the foundation for the advancement of theoretical physics and mathematics "Basis" 17-14-455-1.