Probing the proton structure with associated vector boson and heavy flavor jet production at the LHC

We consider the production of $Z$ bosons associated with heavy (charm and beauty) jets at the LHC energies using two scenarios based on the transverse momentum dependent (TMD) parton densities in a proton. The first of them employs the Catani-Ciafaloni-Fiorani-Marchesini gluon evolution and is implemented in the Monte-Carlo event generator PEGASUS. Here, the heavy quarks are always produced in the hard partonic scattering. The second scheme is based on the parton branching approach, currently implemented into the Monte-Carlo event generator CASCADE. In this scenario, the $Z$ + jets sample is generated and then events containing the heavy flavor jet in a final state are selected. We compare the predictions obtained within these two TMD-based approaches to each other, investigate their sensitivity to the TMD gluon densities in a proton and estimate the effects coming from parton showers and double parton scattering mechanism. Additionally, we compare our predictions with the results of traditional (collinear) pQCD calculations performed at NLO accuracy. It is shown that the TMD-based results agree with the LHC experimental data collected at $\sqrt s = 8$ and $13$ TeV. We discuss the sensitivity of observables to the quark distributions in a proton and present predictions to search for the intrinsic charm signal in forthcoming analyses of the LHC experimental data.

1 Introduction fitted to experimental data, so that the relevant theoretical predictions, where the parton shower effects are already taken into account, can be obtained with no further free parameters -that is in contrast to the usual parton shower event generators 4 .
In the present paper we improve the early calculations [20,21] of associated Z boson and heavy flavor jet production by including into the consideration the effects of parton showers in the initial and final states and extend them to the latest LHC data on Z + c-jet production collected at √ s = 8 and 13 TeV [2,4,6]. The predictions, based mainly on the CCFM gluon dynamics in a proton, will be compared with the results obtained in the PB scenario, implemented in the Monte-Carlo event generator cascade [29]. Such comparison between the calculations performed within these two approaches could be a general consistency check for the k T -factorization phenomenology. At this point, our study is complementary to recent investigations [21,30]. Special interest is related to the comparison of the TMD-based predictions and traditional (collinear) pQCD ones calculated by taking into account higher-order terms. We consider predictions from the standard MadGraph5 amc@nlo tool [31]. Another goal is connected with studying the heavy quark density functions in a proton, which is particularly interesting for the analysis of hard processes at LHC energies. We investigate the influence of IC contributions on various kinematical distributions in Z + c-jet production (and, of course, on the recently measured σ(Z + c)/σ(Z + b) relative production rates). We describe new observables which are sensitive to the IC content of a proton. In this sense we continue the line of our previous studies [32,33]. Finally, we investigate the role of an additional mechanism of Z + c production, double parton scattering (DPS), which is widely discussed in the literature at present (see, for example, [34][35][36][37][38][39][40][41][42] and references therein).
The outline of our paper is following. In Section 2 we briefly describe our theoretical input. The numerical results and discussion are presented in Section 3. Our conclusions are summarized in Section 4.

Theoretical framework
In the present paper to calculate the total and differential cross sections of associated Z boson and heavy flavor jet production at LHC conditions we apply two schemes based on the k T -factorization formalism, which can be considered as a convenient alternative to higher-order DGLAP-based calculations. The first scheme was proposed in [20] and relies mainly on the O(αα 2 s ) off-shell (depending on the transverse momenta of initial particles) gluon-gluon fusion subprocess: which gives the leading contribution to the production cross section in the small x region, where the gluon density dominates over the quark distributions. An essential point here is using the CCFM evolution equation to describe the QCD evolution of the transverse momentum dependent (TMD) gluon density in a proton (see [27]). This equation smoothly interpolates between the small-x BFKL gluon dynamics and high-x DGLAP one, thus providing us with a suitable tool for the phenomenological study. In addition to that, we take into account several subleading subprocesses involving quarks in the initial stateflavor excitation subprocess quark-antiquark annihilation subprocess and quark-gluon scattering which could play a role at large transverse momenta (or, respectively, at large x) where quarks are less suppressed or can even dominate over the gluon density. The last subprocess is taken into account since it provides additional heavy quarks, despite they are obviously suppressed in strong coupling α s . Thus, taking into account the subprocesses (2) -(4) extends the predictions to the whole kinematic range. Note that one has at least one heavy quark Q in the final state already at the amplitude level. The gauge-invariant off-shell amplitude for subprocess (1) was calculated earlier [43,44], where all details are explained. In contrast with the off-shell gluon-gluon fusion, the contributions from quark-involved subprocesses (2) -(4) are taken into account using the DGLAP-based factorization scheme, which provides better theoretical grounds in the region of large x. The evaluation of the corresponding production amplitudes is straightforward and needs no explanation. We only note that the subsequent decay Z → l + l − (including the Z/γ * interference effects) is incorporated already at the production step at the amplitude level in order to fully reproduce the experimental setup. To calculate the contribution from the off-shell gluon-gluon fusion subprocess (1) we used two latest sets 5 of CCFM-evolved TMD gluon densities in a proton, namely, JH'2013 set 1 and set 2 [46]. Their input parameters have been derived from a description of high precision HERA data on proton structure functions F 2 (x, Q 2 ) and/or F c 2 (x, Q 2 ). For quark-induced subprocesses (2) -(4) we have applied the standard CT14 (NNLO) set [47].
The scheme [20] represents a combination of two techniques with each of them being used at the kinematic conditions where it is best suitable. This scheme is implemented into the Monte-Carlo event generator pegasus [48], which has been used in the numerical calculations below. Additionally, we simulate here the effects of parton showers in the initial and final states using the pythia8 [49], thus improving the previous consideration 6 [20,21]. The resulting partons are then processed with fastjet [50] to reconstruct jets in anti-k T algorithm with radia R jet corresponding to the experimental setup. As the heavy quark jet we take the jet, which passes kinematical cuts of the experiment and in which a heavy quark is situated closest to the jet axis in ∆R = ∆η 2 + ∆φ 2 , where ∆η and ∆φ are the corresponding differences in pseudo-rapidity and azimuthal angle. In order to avoid the double counting the during parton shower simulation, we keep the subprocess (4) at parton level calculations only.
We compare our results with a more rigorous scheme based on the Parton Branching (PB) approach [28,29], which provides a solution of the DGLAP equations for conventional and TMD quark and gluon distributions in a proton. The splitting kinematics at each branching vertex is described by the DGLAP equations. Instead of the usual DGLAP ordering in virtuality, angular ordering condition for parton emissions is applied. One of the advantages of this approach is that the PB TMDs can be combined with standard (onshell) production amplitudes, which can be calculated at higher orders. Here we use matrix elements calculated with next-to-leading (NLO) order with MadGraph5 aMC@NLO  Table 1: Basic parameters, used for simulations of associated Z + c-jet production. By default experimental cuts for electrons are shown. Cuts for muons are placed in brackets, if differ.
[31] using the HERWIG6 subtraction terms, which are suitable for combination with PB-TMDs. A special procedure is adopted for the transverse momenta of initial partons: a transverse momentum is assigned according to the TMD density, and then the partonparton system is boosted to its center-of-mass frame and rotated in such a way that only the longitudinal and energy components are nonzero. The energy and longitudinal component of the initial momenta are recalculated taking into account the virtual masses [29]. This method keeps the parton-parton invariant mass exactly conserved, while the rapidity of the partonic system is approximately restored. Similar to the CCFM scenario, the PB TMD parton densities can be obtained via fitting to precise DIS data. Two sets, which differ from each other by a choice of the scale in QCD coupling, were obtained in Ref. [51]. In the numerical calculations below we have used the PB-NLO-HERAI+II-2018 set 2. Technically, we generate a Z + jet(s) sample using cascade and then select events which contain the heavy flavor jet(s) in a final state (see also [30]). This is in contrast to the pegasus calculations, where a heavy flavor jet is always presented in the final state, as explained above.
Finally, we turn to the DPS contribution to Z + c production. We apply the factorization formula [34][35][36][37][38][39][40][41][42]: where σ eff is a normalization constant which incorporates all "DPS unknown" into a single phenomenological parameter. A numerical value σ eff 15 mb was obtained, for example, in recent studies [52][53][54][55][56] from fits to Tevatron and LHC data (see also [57]). This will be taken as the default value throughout the paper. The calculation of inclusive Z boson or charm production cross sections is straightforward and needs no special explanations.

Numerical results
Before presenting results of our calculations let us describe our set of parameters. So, following [61], we apply charm and beauty quark masses m c = 1.4 GeV and m b = 4.75 GeV, mass of Z boson m Z = 91.1876 GeV, its total decay width Γ Z = 2.4952 GeV and sin 2 θ W = 0.23122. As it was mentioned above, we kept n f = 4 active (massless) quark flavors in the pegasus calculations, set Λ (4) QCD = 200 MeV and used two-loop QCD coupling according to [46]. The default renormalization scale was taken to be µ 2 R = m 2 Z . The default factorization scale for the off-shell gluon-gluon fusion subprocess was taken as µ 2 F =ŝ + Q 2 T , where Q T is the net transverse momentum of the initial off-shell gluon pair. This choice is dictated mainly by the CCFM evolution algorithm (see [46] for more information). For quark-induced subprocesses (2) and (3) we keep it equal to the renormalization scale.
The PB calculation with CASCADE [62] were calclualted with m c = 1.47 GeV, m b = 4.5 GeV, α s (m 2 Z ) = 0.118 and µ R = µ F = 1 2 i m 2 i + p 2 t,i , where the sum runs over all particles and parton in the matrix element. The hard process calculations are performed at NLO with MadGraph5 aMC@NLO [31] with herwig6 subtraction terms. The theoretical uncertainties are obtained by varying the scale of the hard process by a factor 2 up and down, provided by MadGraph5 aMC@NLO.
We start from differential cross sections of associated Z + c production at the LHC. Results of our calculations are presented in Figs. 1 -3 in comparison with the CMS data [5,6] taken at √ s = 8 and 13 TeV. The kinematical cuts and jet reconstructing parameters were taken the same as in corresponding experimental analyses (we summarized them in Table 1). The shaded bands represent the theoretical uncertainties of our calculations. To estimate the latter in the pegasus simulation we have used auxiliary "+" and "−" TMD gluon densities in a proton instead of default ones when calculating the off-shell gluon-gluon fusion subprocess (1). These two sets refer to the varied hard scales in the strong coupling α s in the off-shell amplitude: "+" stands for 2µ R , while "−" refers to µ R /2. This was done to preserve the intrinsic consistency of CCFM-based calculations (see [46] for more information). For the quark-induced subprocesses (2) -(4) we just vary the hard scales around its default value between halved and doubled magnitude, as it usually done. We find that the measured Z + c-jet production cross sections are reasonably well reproduced by the pegasus calculations (within the theoretical and experimental uncertainties), although some underestimation of the CMS data taken at √ s = 8 TeV is observed at low p T (c) and large p T (Z). A similar description of the 8 TeV data is achieved in the traditional (collinear) NLO pQCD evaluations, as one can see in Fig. 1. At the same time it is worth pointing out that the two analyses [5,6] used different techniques for the experimental charm jet identification. Namely, in the 8 TeV measurement [5] several methods were utilized for charm identification, including those based on the presence of a muon in the jet or the reconstruction of D ± or D * (2010) ± meson exclusive decays. The 13 TeV measurement benefited however from dedicated machine learning methods developed for identification of charm jets [6]. Interestingly, unlike the pegasus predictions, the cascade results tend to overestimate the 8 TeV data. The predictions of both TMD-based approaches as well as NLO pQCD ones are close to each other at √ s = 13 TeV. The calculated contribution from the DPS production mechanism is small for both considered energies and can play a role at low transverse momenta only.
We find that the pegasus predictions substantially depend on the TMD gluon used, as one can see in Figs. 1 and 2. This can be explained by the fact that the off-shell gluon fusion subprocess (1) plays an essential role in the considered kinematical region. Our calculations show that the CMS data taken at √ s = 8 TeV are better described by JH'2013 set 1 gluon density (except for the first bin at 25 < p T < 30 GeV). Moreover, these predictions practically coincide with the corresponding results of the NLO pQCD calculations. The quark-induced contributions (2) -(4) become important at high transverse momenta, where the typical x values are large, that supports using DGLAP dynamics for these subprocesses. Of course, these subprocesses should be taken into account to describe the data in the whole p T range.
To investigate the effects originating from initial and/or final state parton showers in the scheme implemented into the pegasus tool, we show separately the results obtained at the parton level, that corresponds to the previous calculations [20,21]. We find that   simulation of parton showers leads to some decrease of the calculated cross sections. However, the estimated effect is almost negligible and lies mostly within the bands of theoretical uncertainties.
Concerning the relative σ(Z + c)/σ(Z + b) production rate, we find that the pegasus tends to underestimate recent CMS data, whereas the cascade tool gives better description of the latter. In fact, there is only some underestimation of this ratio at the large transverse momenta of Z boson, see Fig. 3. The observed difference between the pegasus and cascade predictions can be explained by different treatment in the two approaches: in pegasus one always has a heavy flavor jet in the final state, while cascade operates with a sample containing Z+any jets, from which only events having heavy flavor jets after showering are considered [30]. However, the two methods are both compatible with the data within ∼ 2σ. Now we turn to the next point of our study connected with the investigation of heavy quark densities in a proton. In fact, the production of vector bosons accompanied by heavy flavor jets in pp collisions at the LHC can be considered as an additional tool to study the quark and gluon densities in a proton. As it is shown in [12,20], the sensitivity of p t spectra of prompt photons, Z-bosons and c-jets produced in pp → γ/Z + c + X processes at LHC energies to different proton PDFs without the inclusion of the IC component is very small, it is about a few percents. It would be very interesting to study the similar sensitivity to PDFs, which include the IC contribution. To investigate these IC effects in more detail, we repeat the calculations of associated Z + c-jet production cross sections using the CT14 (NNLO) parton densities adopted for BHPS1 (corresponding to the IC probability w max IC = 1%) and BHPS2 (with w max IC = 3.5%) scenarios. Results of our calculations, performed with Monte-Carlo generator pegasus, are shown in Figs. 1 and 2. We find that the IC component is almost undetectable in the kinematical conditions of the CMS experiments even at high transverse momenta. Moreover, even being estimated within the BHPS2 scenario, IC signal lies within the bands of scale uncertainties of our calculations. This agrees with the earlier results [8,12,20,32,33,70], where it was shown that the IC signal can be sizable in the forward rapidity region, |y| ≥ 1.5. The IC effect could be more visible, especially at large transverse momenta (about of 100 GeV and higher) in the relative production rate σ(Z + c)/σ(Z + b) since most of theoretical uncertainties cancels out in this ratio (see Fig. 4).
However, it is known that hadronization effects can result in a significant decreasing of the Z+heavy flavor jet production cross sections at least in some kinematical regions [30]. We check the effect by applying the hadronization effects in our pegasus calculations for the σ(Z + c)/σ(Z + b) cross section ratios measured by CMS [7]. Our results (Fig. 4) show that the hadronization corrections for Z + c production are essentially smaller than for Z + b production. This results in the increasing of the cross section ratio, especially at large transverse momenta of the heavy quark jet (thus leading to better agreement with the data). This agrees with results of [30], obtained with cascade. So the IC effects can be in fact hidden by the hadronization effects, at least for p T (Q) differential cross sections.
The following simple argument can be also useful for further IC studies. Assuming that the IC distribution peaks at x ∼ 0.5, one can expect that the distribution in Feynman variable x F = 2p z / √ s [12,74] would generally follow the initial quark density thus giving rise to an enhancement of the cross sections at large x F values even in a specific kinematical region. Based on this point, we have calculated the cross section of Z + c production as a function of x F using the pegasus tool. The results of our evaluations are shown in Fig. 5. The kinematical cuts applied are given in the last column of Table 1. Note that here we limit ourselves to x F < 0.6 to control statistical uncertainties. We find that, even for the quite conservative IC fraction, the predictions of BHPS1 scenario starts to lie over the uncertainty band of a null hypothesis at x F 0.15. It can illustrate qualitatively the kinematical region, where the IC signal can be visible.

Conclusion
We have considered the production of Z bosons associated with heavy (charm and beauty) jets at the LHC energies using two TMD-based scenarios. The first approach employs the CCFM gluon evolution and has been implemented in the Monte-Carlo event generator pegasus. In this scenario, the heavy quarks are always produced in the hard partonic scattering subprocesses (1) -(4). The second scheme is based on the PB approach implemented into the Monte-Carlo event generator cascade. The traditional NLO pQCD calculations were done also using the standard MadGraph5 amc@nlo tool.
The main goal of this paper is to check the sensitivity of our results to inputs used by calculations, namely: two different TMD gluon distributions (JH'2013 set 1 and JH'2013 set 2), the contribution of the parton showers, the contribution of double parton scattering, different schemes of the QCD calculation, the QCD scale uncertainty, different sets of the conventional PDFs including also the intrinsic charm contribution. We find that there is a sensitivity of transverse momentum distributions of the Z-boson and c-jet to different TMD gluon densities in a proton. In particular, the JH'2013 set 1 gluon leads to a small increase of the p T spectra of both Z boson and c-jet about a few percents in the whole rapidity region. The p T spectra of c-jet or Z boson after inclusion of PS, DPS and IC are changed also by about a few percents at |y| ≤ 2.4. The sensitivity of all our results to the QCD scale uncertainty is about 10 percent in the whole rapidity range. We show also that the contribution from the double parton scattering mechanism is rather small and can play a role at low transverse momenta only. It has been shown that the IC contribution to the ratio σ(Z + c)/σ(Z + b) as a function of heavy flavor jet transverse momentum integrated over the rapidity |y| ≤ 2.4 can be hidden by the hadronization effects. We have illustrated qualitatively that the IC signal can be visible in the x F distribution of c-jet at x F > 0.1, which roughly correponds to large values of p T (Q) and the rapidity range |y| > 1. 5. It has been found that both considered TMD approaches provide a more or less consistent description of recent experimental data on the Z boson and c-jet transverse momentum distributions. This can be seen from a direct comparison between the pegasus and cascade predictions and CMS data collected at √ s = 8 and 13 TeV. Similar agreement with the data is achieved also with MadGraph5 amc@nlo. However, the Monte-Carlo generator cascade provides a better description of the relative σ(Z + c)/σ(Z + b) production rate, that is connected with the different jet production mechanisms implemented into the cascade and pegasus.