Associated non-prompt $J/\psi + \mu$ and $J/\psi + J/\psi$ production at LHC as a test for TMD gluon density

We consider the associated production of $J/\psi$ mesons and muons originating from the $b$-flavored hadron decays and non-prompt double $J/\psi$ production at the LHC using the $k_T$-factorization approach. To describe the inclusive $b$-hadron decays into the different charmonium states we apply fragmentation approach and adopt fragmentation functions based on the non-relativistic QCD factorization. The transverse momentum dependent (TMD) gluon densities in a proton are determined using the Catani-Ciafaloni-Fiorani-Marchesini equation and Kimber-Martin-Ryskin prescription. We investigate the effects coming from parton showers, estimate the double parton scattering contribution and compare our predictions with the first experimental data taken by the ATLAS and LHCb Collaborations at $\sqrt s = 8$~TeV. These data can serve as an additional test for TMD gluon density function in a proton.

Very recently, the ATLAS Collaboration has presented a new measurement of the production of two b-hadrons at √ s = 8 TeV, where one of these b-hadrons decays into a J/ψ meson (with the subsequent decay of the latter into muon pair) and the other decays into µ + X, resulting in three muons in the final state [1]. The kinematic correlations for pairs of b-hadrons, reconstructed via their inclusive decays into J/ψ mesons, have been also measured by the LHCb Collaboration in the forward rapidity region 2 < y J/ψ < 4.5 [2]. These data provide an additional testing ground for perturbative Quantum Chromodynamics (pQCD) predictions of heavy flavour production, especially at small opening angles between the heavy hadrons, where the relevant theoretical uncertainties are rather large. A number of differential cross sections, including different angular correlations between the final decay muons (or rather J/ψ mesons) was measured for the first time, that stimulated us to perform corresponding calculations in the framework of the k T -factorization approach [3,4] and compare these predictions with the ATLAS and LHCb data [1,2]. The k T -factorization approach is mainly based on the Balitsky-Fadin-Kuraev-Lipatov [5] (BFKL) or Catani-Ciafaloni-Fiorani-Marchesini [6] (CCFM) gluon dynamics at small x and has certain technical advantages in the ease of including higher-order radiative corrections that can be taken into account in the form of transverse momentum dependent (TMD) parton distributions 1 . It has become a widely exploited tool and it is of interest and importance to test it in as many cases as possible. Closely related to this is selection of the TMD parton densities best suited to describe data 2 . These tasks form the major goal of our present article. Additionally, we investigate the influence of parton showers on description of the ATLAS and LHCb data [1,2] and estimate the contributions from the double parton scattering (DPS) mechanism, now widely discussed in the literature. The consideration below continues the line of our previous studies [9,10], where we inspect the evolution details of b-quarks fragmenting and decaying into final state charmonia.
The calculations of non-prompt J/ψ meson and a muon associated production (and non-prompt J/ψ + J/ψ production) involve several main ingredients: the cross sections of b-hadron production, the partial widths of their subsequent decays into the different charmonia states and/or semileptonic decays of b-hadrons. Below we collect the previously tested components of the theory and only briefly recall our main points.
First, to calculate the cross sections of inclusive b-hadron production in pp collisions we apply the k T -factorization approach, mainly based on the O(α 2 s ) off-shell gluon-gluon fusion subprocess: where the four-momenta of all particles are given in the parentheses. The corresponding gauge-invariant off-shell (dependent on transverse momenta of the initial gluons) production amplitude was calculated earlier (see, for example, [11] and references therein). Then, b-flavor production cross section can be obtained as a convolution of the off-shell partonic cross sectionσ * gg (x 1 , x 2 , k 2 1T , k 2 2T , µ 2 ) and the TMD gluon distributions in a proton f g (x, k 2 T , µ 2 ): where k i T being the component of the off-shell gluon momentum k i perpendicular to the beam axis (k 2 i = −k 2 iT = 0, i = 1 or 2), x i is the fraction of longitudinal momentum of the colliding proton and µ 2 is the hard scale. The subsequent fragmentation of the produced b quarks into b-hadrons is described with the Peterson fragmentation function [12] with b = 0.0126. The consistency of this setting was shown in a previous paper [10]. We have tested two families of the TMD gluon distribution functions in a proton, which are widely discussed in the literature and often used in appications. So, we used a numerical solution [13] of the CCFM equation (labeled below as JH'2013 family). The CCFM gluon evolution equation provides a suitable tool since it smoothly interpolates between the small-x BFKL gluon dynamics and high-x DGLAP dynamics. Two sets of the TMD gluon densities were determined from the fits to high precision HERA data on the proton structure functions: JH'2013 set 1, which was determined from a fit to inclusive F 2 (x, Q 2 ) data only, and JH'2013 set 2, which was determined from a fit to both F 2 (x, Q 2 ) and F c 2 (x, Q 2 ). As an alternative choice, we applied the TMD gluon density obtained from the Kimber-Martin-Ryskin [14,15] (KMR) prescription. The KMR approach is a formalism to construct the TMD parton (quark and gluon) densities from well-known conventional ones, developed at leading order [14] (LO) and next-to-leading order [15] (NLO). The key assumption of this approach is that the k T -dependence of the TMD parton distributions enters at the last evolution step, so that the usual DGLAP evolution can be used up to this step. For the input, we used Martin-Stirling-Thorn-Watt set [16] (MSTW'2008) at LO and NLO, respectively. The phenomenological consequences of our different choices for the TMD gluon densities in a proton are discussed below 3 .
Next step of our calculations is connected with the description of the J/ψ production from b-hadron decays ("non-prompt" production). We adopted the so-called fragmentation approach, as it was done earlier [10]. In this approach, the calculated b-hadron cross section has to be convoluted with b → J/ψ + X fragmentation function, which is the longitudinal momentum distribution of the J/ψ meson from b-hadron decay, appropriately boosted along the b-hadron flight direction. The latter have been obtained [18] in the framework of the nonrelativistic QCD (NRQCD) factorization [19,20] using the approach [21] and reasonably agree with the CLEO [22] and BABAR [23] measurements. The formalism [18] was implemented into our calculations without any changes (see [10] for more details). To be precise, we employ the asymptotic expression [18] for the bhadron decay distribution differential in the longitudinal momentum fraction z carried by the produced charmonium state, obtained in the limit |p b | m b , where p b and m b are the momentum and the mass of decaying b-hadron. This approximation is valid within 11% and 5% accuracy for |p b | = 10 GeV and 20 GeV, respectively, that is suitable for our phenomenological study. According to the ATLAS and LHCb experimental setup [1,2], we also took into account feed-down contributions from the excited charmonium states, namely, b → χ cJ + X (with J = 0, 1, 2) and b → ψ(2S) + X decays followed by subsequent radiative decays χ cJ → J/ψ + γ and ψ(2S) → J/ψ + γ using the same approach [18]. We set the branching fractions [24]. Finally, J/ψ → µ + µ − decays were generated according to the phase space with the branching fraction B(J/ψ → µ + µ − ) = 5.961% [24].
To produce muons from b-hadron decays, we simulate their semileptonic decay according to the standard electroweak theory with the branching fraction B(b → µ) = 9.8% [24]. We took into account also muons produced in semileptonic cascade decays (the decay of a c-hadron produced in the decay of a b-hadron) with B(b → c → µ) = 8.02% [24]. Other essential parameters, such as renormalization and factorization scales, masses of produced particles are taken exactly the same as in our previous studies [10].
We close the short description of our calculation steps with DPS contributions, where we apply a simple factorization formula (for details see the reviews [25][26][27] and references therein): where σ eff is a normalization constant which incorporates all "DPS unknowns" into a single phenomenological parameter. A similar expression (with an extra factor of 1/2 due to identity of final state particles) is valid in the case of non-prompt J/ψ + J/ψ production. A numerical value of σ eff 15 mb has been obtained from fits to pp and pp data (see, for example, [28]) and will be taken as the default value throughout the paper. Note that one can easily estimate the relative DPS contribution to the total cross section as σ DPS /σ SPS ∼ σ(bb)/σ eff . This value could be non-zero and DPS could contribute significantly in some kinematical regions. Below we will clarify this point for both considered processes.
We discuss first the associated non-prompt J/ψ + µ production at √ s = 8 TeV. The ATLAS Collaboration has measured the corresponding total and differential cross sections in a restricted part of the phase space (fiducial volume). So, each muon was required to have transverse momentum p T > 6 GeV, the two muons originating from the J/ψ decay must have pseudorapidities |η| < 2.3 and the third muon must have |η| < 2.5 [1]. We implemented the experimental setup used by the ATLAS Collaboration in our numerical program. Several normalized differential cross sections have been presented for the first time, namely: transverse momentum of the three-muon system p T (J/ψ, µ), azimuthal separation between the J/ψ and third muon ∆φ(J/ψ, µ), separation between the J/ψ and third muon in the azimuth-rapidity plane ∆R(J/ψ, µ), separation in rapidity between the J/ψ and third muon ∆y(J/ψ, µ), magnitude of the average rapidity of the J/ψ and third muon y boost (J/ψ, µ), mass of the three-muon system M (J/ψ, µ), ratio of the invariant mass of three-muon system to the transverse momentum of three-muon system m µµµ /p µµµ T and its inverse p µµµ T /m µµµ . In some sense, studying of the normalized differential cross sections could lead to a bit more stringent comparison between data and theory due to reduced experimental (mainly systematic) uncertainties.
We confront our predictions with the available data in Figs. 1 -3. The solid histograms represent our central predictions calculated with fixed renormalization µ R and factorization µ F scales at their default values (see [10] for the detailed description of our input), while the shaded regions correspond to scale uncertainties of our predictions. In the case of CCFM-evolved gluon densities (JH'2013 family), to estimate the latter we used the JH'2013 set 1(2)+ and JH'2013 set 1(2)-sets instead of default JH'2013 set 1(2) distributions. These sets represent a variation of the renormalization scale used in the off-shell production amplitude. The JH'2013 set 1(2)+ set stands for a variation of 2µ R , while set JH'2013 set 1(2)-refects µ R /2 (see [13]). To estimate the scale uncertainties of the KMR predictions, we have varied both renormalization and factorization scales around their default values. As one can see, the calculated cross sections (except y boost (J/ψ, µ) spectrum) strongly depend on the TMD gluon density used. A clear difference in shape between the JH'2013 and KMR predictions is observed for angular correlations ∆R(J/ψ, µ) and ∆φ(J/ψ, µ) (see Fig. 1) and distributions on p T (J/ψ, µ) and M (J/ψ, µ) (see Fig. 2). A better description of all these observables is achieved with the KMR family of gluon distributions. There is only small overestimation of the data at low ∆R(J/ψ, µ) and ∆φ(J/ψ, µ) and in the last bins of p T (J/ψ, µ) and M (J/ψ, µ), although the data are rather close to the estimated uncertainty bands. In contrast, both the JH'2013 gluon densities do not reproduce well the measured shape of angular corre-lations: they underestimate the data at low ∆R(J/ψ, µ) and ∆φ(J/ψ, µ) (especially at p T (J/ψ, µ) > 20 GeV) and tend to overestimate the data at ∆φ(J/ψ, µ) ∼ π (see Fig. 1). This is in agreement with the general trend observed in the bb di-jet production [11]. The measured M (J/ψ, µ) distribution is not reproduced with the CCFM-evolved gluon densities. While the y boost (J/ψ, µ) distribution is reasonably well described within the theoretical and experimental uncertainties, the ∆y(J/ψ, µ) spectrum is somewhat poorly described by all gluon densities under consideration with the tendency to fall away at high ∆y(J/ψ, µ). The predictions, obtained with the KMR gluon density, calculated with the LO and NLO accuracy, are close to each other (and even coincide practically within the uncertainties). The estimated DPS contributions are found to be small in the considered kinematic region. Of course, some reasonable variations in σ eff 15 ± 5 mb would affect DPS predictions, though without changing our conclusion.
The observed difference between the KMR and JH'2013 predictions for transverse momentum p T (J/ψ, µ) and invariant mass M (J/ψ, µ) distributions leads to a noticable difference in the m µµµ /p µµµ T and its inverse p µµµ T /m µµµ spectra, as it is shown in Fig. 3. Moreover, this difference becomes even more clearly pronounced: the ratio of the KMR and JH'2013 predictions reaches ∼ 2.5 − 3 at p µµµ T /m µµµ > 2. Therefore, such observables are particularly sensitive to the non-collinear gluon evolution dynamics and, in addition to well-known properties of angular correlations between the momenta of the produced particles, could be very promising to constrain the TMD gluon densities in a proton. Now we turn to the non-prompt J/ψ + J/ψ production. The LHCb Collaboration has presented [2] the normalized cross sections measured as functions of several variables, namely, the transverse momentum, rapidity and invariant mass of the J/ψ pair, the difference in the azimuthal angle between the momentum directions of two J/ψ mesons, the difference in the rapidity and pseudorapidity between them and the assymetry A T between the transverse momenta of produced mesons: These data were obtained in the forward rapidity region 2 < y J/ψ < 4.5 for different requirements on the minimum transverse momentum of the J/ψ mesons. The results of our calculations are shown in Figs. 4 -10 in comparison with the LHCb data [2]. One can see that, in general, all the TMD gluon densities under consideration describe the data reasonably well for all distributions within the uncertainties. However, both the CCFM-evolved gluons tend to overestimate the measured transverse momentum spectra of the J/ψ pair at low p T (J/ψ, J/ψ) and underestimate the invariant mass distributions near the threshold, where M (J/ψ, J/ψ) ≤ 10 GeV. The KMR calculations agree well with the LHCb data for these observables. We see again that a clear difference between the predictions is observed in the angular correlations (see Fig. 9). Similar to the J/ψ + µ production, both the JH'2013 gluon densities do not reproduce the measured shape of ∆φ(J/ψ, J/ψ) distributions: they underestimate the LHCb data at low ∆φ(J/ψ, J/ψ) and overestimate the data at ∆φ(J/ψ, J/ψ) ∼ π, although at large transverse momenta overall agreement becomes better. In contrast, the both KMR gluons provide good description of these angular correlations.
The DPS contributions are also found to be small in the considered kinematic region. Note that it was discussed earlier [29] that a special correction factor should be included into DPS calculations at LHCb conditions to take into account limited partonic phase space. We have found that such a factor F ∼ (1 − x 1 − x 2 ) 2 results in ∼ 15 − 20% smaller DPS cross sections (not shown in Figs. 1 -10), that, in any case, does not change our conclusions. For the ATLAS kinematical region the influence of this correction factor is negligible.
As a last point of our study, we would like to note that the angular correlations in bb production (and, therefore, ∆φ(J/ψ, µ), ∆R(J/ψ, µ) and ∆φ(J/ψ, J/ψ) distributions considered above) are known to be sensitive to the inclusion of (initial and final state) parton radiation that can either be simulated as parton showers or taken as an additional higher-order process 4 . For the latter, we have estimated the partial contribution coming from the O(α 3 s ) quark-gluon scattering subprocess qg → qbb within the conventional (collinear) QCD factorization using the madgraph tool [30]. We found it to be negligible everywhere (not shown in the Figs. 1 -10). Additionally, we investigated the influence of parton showers on the description of the ATLAS and LHCb data [1,2]. For these studies we used a parton shower algorithm implemented in the Monte-Carlo event generator cascade [31]. The results of our calculations for ∆φ(J/ψ, µ), ∆R(J/ψ, µ) and ∆φ(J/ψ, J/ψ) distributions are shown in Figs. 11 and 12, where the JH'2013 set 2 gluon density is applied as an example. We observed only a small effect at low ∆φ(J/ψ, µ) and ∆φ(J/ψ, J/ψ) in the kinematical region of ATLAS and LHCb experiments [1,2]. This effect is coming from the final state parton showers only, that could be easily understood since in our calculations the initial state parton emissions are already determined from the TMD gluon density and, therefore, do not influence the k T of the gluons. However, parton showers could be very important at some kinematical regions (see, for example, [32,33]).
To conclude, we applied the k T -factorization approach to investigate the associated production of J/ψ mesons and muons originating from the b-hadron decays and nonprompt J/ψ + J/ψ production at the LHC. Our calculations were inspired by the recent experimental data presented by the ATLAS and LHCb Collaborations at √ s = 8 TeV, where a number of differential cross sections were measured for the first time. To describe the inclusive b-hadron decays into the J/ψ mesons we used fragmentation approach with corresponding fragmentation functions calculated within the NRQCD. Following the experimental setup, we took into account both direct J/ψ production mechanism and feed-down contributions from the radiative decays of excited charmonium states. The semileptonic b-hadron decays were generated according to the standard electroweak theory. Numerically, we have tested two families of the TMD gluon densities in a proton, namely, the CCFM-evolved gluon distributions (JH'2013 sets) and the KMR ones, calculated with the LO and NLO accuracy. Theoretical uncertainties and effects arising from parton showers, higher-order pQCD corrections and double parton scattering mechanism were estimated. Quite satisfactory agreement between the predictions and the data can be regarded as another voice in support of the chosen TMD gluon parametrizations. The tri-muon and J/ψ + J/ψ measurements at the LHC did and will continue to play their role in providing the useful and necessary experimental constraints.        no PS with PS ATLAS Figure 11: Influence of the parton shower effects on the angular correlations in the associated non-prompt J/ψ + µ production at √ s = 8 TeV. The JH'2013 set 2 gluon density is applied. The experimental data are from ATLAS [1].  Figure 12: Influence of the parton shower effects on the angular correlations in the nonprompt J/ψ + J/ψ production at √ s = 8 TeV. The JH'2013 set 2 gluon density is applied. Notation of histograms is the same as in Fig. 11. The experimental data are from LHCb [2].