Charm fragmentation and associated $J/\psi + Z/W^\pm$ production at the LHC

We consider the production of electroweak $Z$ or $W^\pm$ bosons associated with $J/\psi$ mesons at the LHC conditions. Our attention is focused on new partonic subprocesses which yet have never been considered in the literature, namely, the charmed or strange quark excitation subprocesses followed by the charmed quark fragmentation $c\to J/\psi+c$. Additionally we take into account the effects of multiple quark and gluon radiation in the initial and final states. We find that the contributions from the new mechanisms are important and significantly reduce the gap between the theoretical and experimental results on the $J/\psi + Z$ and $J/\psi + W^\pm$ production cross sections.


Motivation
Hadronic production of J/ψ mesons in association with electroweak gauge bosons (Z or W ± ) are interesting processes [1,2]. They involve both strong and weak interactions and may serve as a complex test of perturbative QCD, electroweak theory and parton evolution dynamics. Moreover, they provide a unique laboratory to investigate the charmonia production mechanisms predicted by the non-relativistic QCD (NRQCD) factorization [3,4]. The latter is a rigorous framework for the description of heavy quarkonia production and/or decays and implies a separation of perturbatively calculated short distance cross sections for the production of a heavy quark pair in an intermediate Fock state 2S+1 L (a) J with spin S, orbital angular momentum L, total angular momentum J and color representation a from its subsequent non-perturbative transition into a physical quarkonium via soft gluon radiation. However, the charmonia production at high energies is not fully understood at the moment despite many efforts made in last decades. In fact, NRQCD has a long-standing challenge in the J/ψ and ψ polarization and provides inadequate description of the η c production data 1 (see, for example, [6][7][8] for more information). The associated production of J/ψ mesons and gauge bosons at the LHC is a suitable process for further studying NRQCD.
A complete next-to-leading order (NLO) predictions for prompt J/ψ + Z/W ± production within the NRQCD are available [9][10][11]. It was shown that the differential cross sections at the leading order (LO) are significantly enhanced by the NLO corrections. There are only color octet (CO) contributions to the J/ψ +W ± production at both LO and NLO level [9,10], while color singlet (CS) terms contribute at higher orders [9]. Taken together with the corresponding predictions from the double parton scattering (DPS) mechanism, the NLO NRQCD expectations for associated J/ψ + Z production cross sections are lower than the data [1] by a factor of 2 to 5 (depending on the J/ψ transverse momentum). The difference between the theoretical predictions and the measured J/ψ + W ± cross sections is typically even larger [2]. Of course, these discrepancies need an explanation. So, further theoretical studies are still an urgent and important task.
In the present paper, we draw attention to a new contributions to prompt J/ψ+Z/W ± production cross section, namely, the flavor excitation subprocesses charm quark excitation for Z bosons (or strange quark excitation for W ± bosons) followed by subsequent charm fragmentation, c → J/ψ + c. These contributions have been overlooked in the literature and yet have never been considered. Of course, this flavor excitation is of formally non-leading order in α s in the pQCD expansion, but it could play a role because of its different momentum dependence (in comparison with what was previously considered [9][10][11]). The main goal of our study is to clarify this point. In addition, we calculate contributions coming from multiple gluon radiation in the initial state and final states through CO gluon fragmentation (via g → cc[ 3 S 1 ] → J/ψ channel). This mechanism was found to be important [12] for double J/ψ production in the kinematic region covered by the CMS and ATLAS measurements. We see certain interest in examining this kind of contributions for other processes, such as vector boson production in association with J/ψ mesons.
The outline of the paper is the following. In Section 2 we describe the basic steps of our calculations. In Section 3 we present the numerical results and discussion. Our conclusions are summarised in Section 4. Figure 1: Example of Feynman diagram taken into account in the NRQCD calculations [9][10][11] (left panel) and diagram of charm excitation followed by the c-quark fragmentation to J/ψ (right panel).

The model
As it was already mentioned above, the new subprocesses can be described as the sea (charm or strangeness) excitation followed by the c-quark fragmentation: and similar subprocesses forc ands antiquarks. There is no double counting with the previous estimations [9][10][11] based on the quark-antiquark annihilation and gluon-gluon fusion subprocesses since the subprocesses (1) have rather different final state containing c-quark (as it is clearly seen in Fig. 1). Despite non-leading in α s , the subprocesses (1) have some advantage in kinematics. The large mass m of the emitted boson in (2) suppresses two quark propagators and leads to the dependence σ ∼ 1/m 8 , see Fig. 1 (left diagram). In contrast, the behavior of quark excitation subprocesses (1) is σ ∼ 1/m 4 . This compensates the sparseness of the charmed or strange sea and the presence of extra coupling α s . In the calculations shown below, we exploit the idea that the sea quarks all appear from a perturbative chain as a result of the QCD evolution of gluon densities in a proton. Then, appending an explicit gluon splitting vertex to (1), we come to the gluon-gluon fusion subprocesses, see Fig. 2 (a) and (b): The W + +c + s production subprocesses can be obtained by charge conjugation. The difference between (1) and (3) is that the sea quark density in (1) may contain a nonperturbative component. Yet another class of new contributions is represented by the initial and final state gluon radiation that accompanies the production of electroweak bosons, where the radiated gluons then convert into J/ψ mesons via CO channel g → cc[ 3 S 1 ] → J/ψ (see also [12]). Here we take into account the following subprocesses: where q denotes any light quark flavor. The initial and final state parton emission is simulated using the standard parton showering algorithms (see discussion below). To evaluate the subprocesses (3) and (4) we employ the k T -factorization approach [13,14], which can be considered as a convenient alternative to explicit fixed-order perturbative QCD calculations at high energies. We see certain technical advantages in the fact that, even with  [18] dynamics, thus providing us with the suitable tool for our phenomenological study. The gauge-invariant off-shell production amplitudes for subprocesses (3) and (4) have been calculated in [19,20] and all of the technical details are explained there 3 . To reconstruct the CCFM evolution ladder, we generate a Les Houches Event file using the Monte-Carlo event generator pegasus [21] and then process the file with the TMD parton shower routine implemented into the Monte-Carlo event generator cascade [22]. In this way we obtain full information about gluon emissions in the initial state.
For completeness, we also take into account several subprocesses involving quarks in the initial state shown in Fig. 2 where we keep only valence quarks in the first two cases to avoid double counting with off-shell gluon-gluon fusion subprocesses (4), see Fig. 2 (c) and (d). The last subprocesses, shown in Fig. 2 (e) -(h), are taken into account since they provide additional charmed quarks, despite they are obviously suppress in α s . The quark-induced contributions (5) can play a role at essentially large transverse momenta (or, respectively, at large x which is needed to produce high p T events) where the quarks are less suppressed or can even dominate over the gluon density. Here we find it reasonable to rely upon the conventional (collinear) DGLAP-based factorization, which provides better theoretical grounds in the region of large x. So, our scheme represents a combination of two techniques with each of 2 A detailed description of this approach can be found, for example, in review [15]. 3 These amplitudes are implemented in the Monte-Carlo event generator pegasus [21]. them being used at the kinematic conditions where it is best suitable 4 . The evaluation of the relevant hard scattering amplitudes and production cross sections is straightforward and needs no explanation. To simulate the parton emissions in initial and final states, we employ the parton shower routine implemented into the Monte-Carlo event generator pythia8 [26]. In this way, we collect all the radiated gluons and charmed quarks.
The parton-level amplitudes (3) -(5) provide the starting point for the charmed quark and the gluon fragmentation. The reliability of the fragmentation approach can be motivated by the high transverse momenta typical for the recent ATLAS data [1,2] taken at √ s = 8 TeV (see, for example, [27,28]). The fragmentation function D H a (z, µ 2 ) describes the transition of a parton a into a charmonium state H and can be expressed as a collection of contributions from the different intermediate states: where n labels the intermediate CS or CO state of the charmed quark pair, µ 2 is the fragmentation scale and O H [n] are the corresponding long-distance matrix elements (LDMEs) [3,4]. Then, the eventual cross section reads where p (c) , p (g) and p are the momenta of the c-quark, gluon and outgoing charmonium state H, respectively. Only the c → cc[ 3 S 1 ]+c and g → cc[ 3 S 1 ] transitions give sizeable contributions to the S-wave charmonia (J/ψ, ψ ). For P -wave mesons χ cJ with J = 0, 1 or 2 one has to take into account the leading terms c → cc[ 3 P In our numerical calculations we take into account the feeddown contributions from χ c1 , χ c2 and ψ decays and neglected the χ c0 decays because of low branching fraction of the latter [29].
The scale dependence of the fragmentation functions is driven by the multiple gluon radiation. We take the initial conditions in the form (see, for example, [30]) where the starting scale is µ 2 0 = m 2 ψ . The effects of final state multiple radiation can be described by the LO DGLAP evolution equation where P ab are the standard LO DGLAP splitting functions. According to the NRQCD approximation, we set the charm mass to m c = m H /2 and then solve the DGLAP equations (12) numerically. Concerning other parameters involved into our calculations, we have used the TMD gluon densities in a proton obtained from the numerical solution of CCFM evolution equation, namely, JH'2013 set 1 and set 2 gluons [31]. The input parameters of JH'2013 set 1 gluon distribution have been fitted to the precise HERA measurement of the proton structure function F 2 (x, Q 2 ), whereas the input parameters of JH'2013 set 2 gluon were fitted to the both structure functions F 2 (x, Q 2 ) and F c 2 (x, Q 2 ). According to [31], we use the two-loop formula for the QCD coupling α s with n f = 4 active quark flavors at Λ 1 ] = 0.7038 GeV 3 (see, for example, [3,4] and references therein). The CS wave functions at the origin for χ cJ mesons and CO LDMEs for charmonia family have been determined [33,34]:

Numerical results
In this section, we present the results of our calculations and perform a comparison with the recent ATLAS data [1,2]. The ATLAS Collaboration has measured the differential cross sections of associated Z/W ± + J/ψ production as a function of J/ψ transverse momentum in a restricted part of the phase space (fiducial volume) at √ s = 8 TeV. In the case of Z + J/ψ production, the J/ψ meson was required to have transverse momentum p J/ψ T > 8.5 GeV and rapidity |y J/ψ | < 2.1, while the leading and subleading muons originated from subsequent Z boson decays must have pseudorapidities |η l | < 2.5 and transverse momenta p l T > 25 GeV and p l T > 15 GeV, respectively. Their invariant mass m ll is required to be |m ll − m Z | < 10 GeV, where m Z is the Z boson mass. For W ± + J/ψ production, the following cuts are applied: p J/ψ T > 8.5 GeV, |y J/ψ | < 2.1, p l T > 25 GeV and |η l | < 2.4 for muon and p ν T > 20 GeV for neutrino originated from the W ± boson decays. The transverse mass of W ± boson defined as is required to be m T (W ± ) > 40 GeV, where φ l and φ ν are the azimuthal angles of the decay muon and neutrino. We have implemented the experimental setup used by the ATLAS Collaboration in our calculations. Our numerical results are shown in Fig. 3, where we have used JH'2013 set 1 gluon density as the default choice. We find that the contributions from subprocesses (3) -(5) with their subsequent parton fragmentation into J/ψ mesons are remarkably important, especially at large transverse momenta. In fact, at p J/ψ T ≥ 20 − 30 GeV it gives approximately the same contribution as that coming from the NLO NRQCD estimations summed with the corresponding contribution from the DPS production mechanism (we took the latter from the ATLAS papers [1,2]). The contributions (3) -(5) are large for W ± + J/ψ production (where only the CO contributions are presented in the NLO NRQCD predictions) almost in the whole range of p J/ψ T . One can see that summing the contributions from subprocesses (3) -(5) and NLO NRQCD predictions (and DPS terms, of course) allows us to significantly reduce the discrepancy between the theoretical expectations and experimental data. Moreover, the upper edge of our estimated uncertainty band shown in Fig. 3 is rather close to the ATLAS data for both processes under consideration.
As usual, to evaluate the uncertainties we have varied the renormalization and fragmentation scales around their default values by a factor of 2. We only note that we replaced the JH'2013 set 1 gluon density by the JH'2013 set 1+ or JH'2013 set 1− ones when calculating the uncertainties connected with the variation of renormalization scale in (3) and (4). This was done to preserve the intrinsic consistency of the calculation, that is, to observe the correspondence between the TMD gluon set and scale used in the CCFM evolution (see [31] for more information). The dominant uncertainties are found to come from the scales of gluon densities, while the ones coming from fragmentation functions are almost negligible compared to the first.
To investigate the sensitivity of our results to the choice of TMD gluon density we repeat the calculations with JH'2013 set 2 distribution. We find that the obtained predictions are rather close to the ones obtained with the default JH'2013 set 1 gluon (or even coincide with them within the uncertainties). The feeddown contributions from ψ and χ cJ decays also play a significant role. Their contribution is about of 30% of the estimated direct contribution in a wide p J/ψ T range, as one can see in Fig. 3. Now let us discuss the role of multiple parton radiation. As it was mentioned above, we considered here two qualitatively different sources of parton fragmentation into the J/ψ mesons, namely, fragmentation of charmed quarks, originated in the hard interaction and fragmentation of gluons, originated as a result of initial QCD evolution of parton densities. Note that the charmed quarks emitted in the initial state according to pythia parton showering algorithm were attributed to the charm quark fragmentation to keep the charm/gluon separation. The effects of parton showers in the final states are taken into account in the form of DGLAP-evolved fragmentation functions. The contributions from charmed quarks and gluon fragmentation are shown in Fig. 4 separately. In the case of gluon-gluon fusion subprocesses, the main contribution at the small transverse momenta p J/ψ T comes from fragmentation of charmed quarks while at the large p J/ψ T the gluon fragmentation starts to dominate. For quark-induced subprocesses (5) small p J/ψ T region is driven mainly by the fragmentation of multiple gluon radiation, whereas the region of high p J/ψ T is governed by the charm quark fragmentation. Nevertheless, as one can clearly see, in both cases the fragmentation of multiple gluon emission noticeably enhances the charm fragmentation and provides a sensible growth of the total and differential cross sections.
In summary, we can conclude that adding the new production mechanisms (3) -(5) to the conventional NLO NRCQD predictions greatly improves the situation for prompt Z/W ± + J/ψ production and significantly decreases the gap between the theoretical estimations and the data. Further on, an appropriate treatment of the initial and final state parton emission and including the feeddown from the χ c and ψ decays enhances the theoretical expectations even more and makes the agreement with the data even better.

Conclusion
We have considered the production of electroweak Z or W ± bosons associated with J/ψ mesons in pp collisions at the LHC at √ s = 8 TeV. We have investigated the role of new partonic subprocesses which yet have never been considered in the literature, namely, the flavor (charm or strangeness) excitation subprocesses followed by the charm fragmentation c → J/ψ +c. In addition, we have taken into account the effects of multiple quark and gluon radiation in the initial and final states. We have demonstrated that the considered new contributions are remarkably important and significantly reduce the gap between the theoretical and experimental results on the J/ψ + Z/W ± production cross sections.