NRQCD Confronts LHCb Data on Quarkonium Production within Jets

TL;DR

This paper compares NRQCD-based theoretical calculations with LHCb data on $J/\psi$ production within jets, showing NRQCD models outperform default PYTHIA simulations in matching the observed $z(J/\psi)$ distribution.

Contribution

It demonstrates that NRQCD calculations, combined with resummation techniques, better describe LHCb jet data than standard PYTHIA simulations, highlighting the importance of nonperturbative effects.

Findings

NRQCD models agree better with LHCb data than PYTHIA.

Resummed FJF formalism provides accurate $z(J/\psi)$ predictions.

High-$p_T$ LDME fits align well with experimental results.

Abstract

We analyze the recent LHCb measurement of the distribution of the fraction of the transverse momentum, $z(J/\psi)$, carried by the $J/\psi$ within a jet. LHCb data is compared to analytic calculations using the fragmenting jet function (FJF) formalism for studying $J/\psi$ in jets. Logarithms in the FJFs are resummed using DGLAP evolution. We also convolve hard QCD partonic cross sections, showered with PYTHIA, with leading order Non Relativistic Quantum Chromodynamics (NRQCD) fragmentation functions and obtain consistent results. Both approaches use Madgraph to calculate the hard process that creates the jet initiating parton. These calculations give reasonable agreement with the $z(J/\psi)$ distribution that was shown to be poorly described by default PYTHIA simulations in the LHCb paper. We compare our predictions for the $J/\psi$ distribution using various extractions of nonperturbative NRQCD long-distance matrix elements (LDMEs) in the literature. NRQCD calculations agree with LHCb data better than default PYTHIA regardless of which fit to the LDMEs is used. LDMEs from fits that focus exclusively on high transverse momentum data from colliders are in good agreement with the LHCb measurement.

Figures

Tables

Table 1. Table 1: LDMEs for NRQCD production mechanisms used in this paper in units of GeV 3 superscript GeV 3 {\rm GeV^{3}} .

	$⟨{𝒪}^{J/\psi }{(}^3{S}_1^{\left[1\right]})⟩$	$⟨{𝒪}^{J/\psi }{(}^3{S}_1^{\left[8\right]})⟩$	$⟨{𝒪}^{J/\psi }{(}^1{S}_0^{\left[8\right]})⟩$	$⟨{𝒪}^{J/\psi }{(}^3{P}_0^{\left[8\right]})⟩/m_c^2$
	$\times$ GeV3	$\times {10}^{-2}$ GeV3	$\times {10}^{-2}$GeV3	$\times {10}^{-2}$GeV3
B & K Butenschoen and Kniehl (2011, 2013)	$1.32\pm 0.20$	$0.224\pm 0.59$	$4.97\pm 0.44$	$-0.72\pm 0.88$
Chao, et al. Chao et al. (2012)	$1.16\pm 0.20$	$0.30\pm 0.12$	$8.9\pm 0.98$	$0.56\pm 0.21$
Bodwin et al. Bodwin et al. (2014)	$1.32\pm 0.20$	$1.1\pm 1.0$	$9.9\pm 2.2$	$0.49\pm 0.44$

Equations

$$d\sigma [pp\to J/\psi X]=\sum _{n}d\sigma [pp\to c\bar{c}(n)X]⟨{\mathcal{O}}^{J/\psi }(n)⟩,$$

$$d\sigma [e^{+}{e}^{-}\to jets,h]=H\otimes S\otimes J(\otimes J)\otimes {\mathcal{G}}^h$$

$$r(i,n)=\frac{d\hat{\sigma }(pp\to i+X){\int }_0^{1}dz{D}_{i\to J/\psi }^n(z)}{d\hat{\sigma }(pp\to c+X){\int }_0^{1}dz{D}_{c\to J/\psi }^{{}^3{S}_1^{[1]}}(z)}.$$

Full text

The production of quarkonium is a challenging test of Quantum Chromodynamics due to the mutiple length scales involved. The LHCb collaboration Aaij et al. (2017) published the first study of $J/\psi$ produced within jets. The distribution of the fraction of the jet’s transverse momentum, $p_{T}$, carried by the $J/\psi$, $z(J/\psi)$, was found to disagree significantly with predictions from the PYTHIA monte carlo Sjostrand et al. (2006); Sjo (2015) using leading order calculations of $J/\psi$ production in the Non-Relativistic Quantum Chromodynamics (NRQCD) factorization formalism Bodwin et al. (1995). This letter is provides improved theoretical calculations of the $z(J/\psi)$ distribution and to discuss the implications of the LHCb results for the NRQCD factorization formalism.

Production of quarkonium in hadron colliders has been the subject of experimental and theoretical studies for decades. The problem is challenging because it involves several disparate scales. These include $p_{T}$, which can be much larger than the mass of the bound state, $\approx 2m_{Q}$, where $m_{Q}$ is the mass of the heavy quark, as well as scales that are much smaller: the relative momenta, $m_{Q}v$ ($v$ is the typical velocity of the heavy quarks in the bound state), the kinetic energy, $m_{Q}v^{2}$, and the nonperturbative scale $\Lambda_{QCD}$.

The most common approach to calculating quarkonium production is the NRQCD factorization formalism Bodwin et al. (1995). In this formalism, the cross section for $J/\psi$ in a $pp$ collision is written as

[TABLE]

where $d\sigma[pp\to c\bar{c}(n)X]$ is the short distance cross section for producing the $c\bar{c}$ pair in a state $n$ with definite color and angular momentum quantum numbers and $\langle{\cal O}^{J/\psi}(n)\rangle$ is a long distance matrix element (LDME) that describes the nonperturbative transition of the $c\bar{c}$ pair in the state $n$ into a final state containing $J/\psi$. $X$ denotes other possible particles in the final state. The quantum numbers $n$ will be denoted ${}^{2S+1}L_{J}^{[i]}$ where the notation for angular momentum is standard and $i=1\,(8)$ for color-singlet (color-octet) states. The short distance cross sections are perturbatively calculable in a power series in $\alpha_{s}$, while the LDMEs are nonperturbative and must be extracted from data. The LDME scale with definite powers of $v$ so the NRQCD factorization formalism organizes the calculation of quarkonium production (and decay) into a systematic double expansion in $\alpha_{s}$ and $v$.

For $J/\psi$ production, the leading matrix element in the $v$ expansion is ${}^{3}S_{1}^{[1]}$ which scales as $v^{3}$. The next most important are the color-octet LDMEs: ${}^{3}S_{1}^{[8]}$, ${}^{1}S_{0}^{[8]}$, and ${}^{3}P_{J}^{[8]}$, which all scale as $v^{7}$. $J/\psi$ production has been measured in a wide variety of experiments, including $e^{+}e^{-}$, $pp$, $p\bar{p}$, $ep$, $\gamma p$, and $\gamma\gamma$ collisions, spanning a wide range of energies. At present, next-to-leadiing order (NLO) QCD calculations are available for the above mentioned color-singlet and color-octet mechanisms for all these initial states. Global fits to the world’s data using these calculations were performed in Refs. Butenschoen and Kniehl (2011, 2013). The resulting LDMEs are shown in the first line in Table 1. The LDMEs are consistent with the expected $v^{4}$ suppression of the color-octet mechanisms. The global fits in Refs. Butenschoen and Kniehl (2011, 2013) are reasonably well described by NLO NRQCD, but there are nagging discrepancies that call into question our understanding of quarkonium production. The most notable discrepancy is the polarization puzzle: if the LDMEs of Refs. Butenschoen and Kniehl (2011, 2013) are used, the produced $J/\psi$ are predicted to be polarized transverse to their momentum at high $p_{T}$, while in fact they are are produced with essentially no polarization. (This is also a problem in $\Upsilon$ production.) Another important discrepancy is the failure of spin symmetry predictions for $\eta_{c}$ production Aaij et al. (2015); Butenschoen et al. (2015). (For possible solutions to the $\eta_{c}$ problem using different extractions of LDMEs see  Han et al. (2015); Zhang et al. (2015); Sun and Zhang (2015).) In light of the failure of NLO QCD to predict the $J/\psi$ polarization, other authors have proposed alternative approaches to fitting the LDMEs. Refs. Chao et al. (2012); Bodwin et al. (2014) have emphasized that NRQCD factorization should be most reliable at the highest values of $p_{T}$ and have performed fits that focus exclusively on high $p_{T}$ $J/\psi$ production in colliders. Ref. Bodwin et al. (2014) also merges NLO calculations with fragmentation contributions in which Altarelli-Parisi evolution is used to resum logs of $p_{T}/m_{J/\psi}$. The LDMEs from the fits of Refs. Chao et al. (2012); Bodwin et al. (2014) are shown in the second and third lines of Table 1, respectively. We will use these three sets of LDMEs in our analysis. There have been other fits to the LDMEs Gong et al. (2013); Bodwin et al. (2016) which include explicit feeddown from $\chi_{cJ}$ states. Since these effects are not included in our calculations, we do not use these LDME extractions in this work.

Recently, Ref. Baumgart et al. (2014) proposed studying the distribution of quarkonia produced within jets as an alternative test of NRQCD in hadron colliders. Cross sections for jets with identified hadrons are given in terms of fragmenting jet functions (FJF) that were first introduced in Ref. Procura and Stewart (2010) and studied further in Refs. Liu (2011); Procura and Waalewijn (2012); Jain et al. (2012, 2011, 2013); Bauer and Mereghetti (2014); Ritzmann and Waalewijn (2014); Kaufmann et al. (2015); Kang et al. (2016, 2017a); Dai et al. (2016, 2017). The FJFs are functions of the jet energy, $E$, and the fraction of energy carried by the identified hadron, $z$. FJFs are calculable as a convolution of the more inclusive fragmentation function with a perturbative matching coefficient evaluated at the jet energy scale, $E_{J}=2E\tan(R/2)$. Ref. Baumgart et al. (2014) showed that the quarkonium FJF can be calculated in terms of NRQCD fragmentation functions Braaten and Fleming (1995); Braaten and Yuan (1993); Braaten et al. (1993) and that the $z$ and $E$ dependence of these cross sections are sensitive to the underlying production mechanisms because NRQCD fragmentations differ for different production mechanisms. For further work see Refs. Bain et al. (2016a, b); Chien et al. (2016).

Ref. Bain et al. (2016a) used the FJF formalism to compute cross sections for jets with $B$ mesons and $J/\psi$ produced within jets in $e^{+}e^{-}$ collisions. For $B$ mesons the paper studied $e^{+}e^{-}\to b\bar{b}$ followed by $b$ quark fragmenting to a jet with a $B$ meson. For $J/\psi$, Ref. Bain et al. (2016a) studied $e^{+}e^{-}\to b\bar{b}g$ followed by gluon fragmentation to a jet with $J/\psi$. In both cases Ref. Bain et al. (2016a) studied the dependence of the cross section on $z$, the fraction of the energy carried by the identified hadrons, and the jet’s angularity, $\tau_{a}$, Berger et al. (2003) a jet substructure variable whose definition can be found in Ref. Bain et al. (2016a). The analytic expression for these cross sections has the schematic form

[TABLE]

where $h$ is either a $B$ meson or $J/\psi$, ${\cal G}^{h}$ is the FJF for the hadron $h$, $J$ is a jet function for the other jets in the event (there is one $J$ for $B$ mesons and two for $J/\psi$), $S$ is the soft function, and $H$ is the hard cross section for $e^{+}e^{-}\to b\bar{b}$ or $e^{+}e^{-}\to b\bar{b}g$. Dependence on all kinematic quantities has been suppressed. Evaluation of each of the quantities appearing in Eq. (1) shows that they all have logarithms of different scales. The renormalization group equation (RGE) for each of the functions in Eq. (1) needs to be solved, and the functions need to be evolved to a common scale so that large logarithms are resummed to all orders in perturbation theory. Details of this formalism for jet cross sections can be found in Ref. Ellis et al. (2010). Analytic calculations in Ref. Bain et al. (2016a) were performed to next-to-leading-log-prime (NLL’) accuracy.111NLL’ means that in addition to NLL accuracy, $O(\alpha_{s})$ contributions to the soft, jet, and FJF are also kept. The calculations were also performed using the PYTHIA monte carlo. For $B$ mesons PYTHIA and NLL’ analytical calculations were in good agreement. However, for jets with $J/\psi$ good agreement was found only in the $\tau_{a}$ distributions. The $z$ distributions predicted by PYTHIA at LO were significantly harder than the $z$ distributions predicted by the NLL’ calculations. This discrepancy between theoretical calculations and PYTHIA is remarkably similar to the discrepancy between data and PYTHIA recently found by the LHCb collaboration. This motivates us to perform calculations similar to those of Ref. Bain et al. (2016a) to obtain a better description of the LHCb data.

To understand the discrepancy between analytical NLL’ calculations of the $z$ distributions and PYTHIA, one must understand how PYTHIA models the production of quarkonium. In PYTHIA the heavy quark-antiquark pair is produced in the short-distance process in either a color-octet or color-singlet state. If it is in a color-singlet state the heavy quark antiquark pair behaves like a color-singlet particle, emits no gluon radiation, and eventually turns into the quarkonium. If the heavy quark-antiquark pair is produced in a color-octet state, PYTHIA treats the pair as a single colored particle that showers with the splitting function $2P_{qq}(z)$. As this splitting function is strongly peaked at $z=1$, the color-octet pair retains most of its momentum after the shower. At the end of the shower the color-octet quark antiquark emits a soft gluon in order to become a color-singlet quarkonium.

The physical picture of quarkonium production implied by the NLL’ analytic calculations in Ref. Bain et al. (2016a) is quite different from PYTHIA. The FJF that controls the $z$ dependence of the cross section is, up to $O(\alpha_{s}(E_{J}))$ corrections, equal to the fragmentation function evaluated at the scale $E_{J}$. The evolution of the fragmentation function from the scale $2m_{c}$ up to the scale $E_{J}$ is governed by Altarelli-Parisi evolution equations. Ref. Bain et al. (2016a) showed that this is equivalent to producing a hard gluon in the short-distance process with virtuality of order $E_{J}$, allowing that gluon to shower until a gluon with virtuality $\sim 2m_{c}$ hadronizes into the $J/\psi$. This can be implemented in PYTHIA by simulating events in which the gluon is produced in the hard process, hadronization is turned off, and allowing the parton to shower down to a scale $\sim 2m_{c}$. After this a gluon $z$ distribution is obtained and convolved manually with a perturbative NRQCD fragmentation function (calculated at lowest order in $\alpha_{s}(2m_{c}))$. This procedure was referred to Gluon Fragmentation Improved PYTHIA (GFIP) in Ref. Bain et al. (2016a), and was shown to give good agreement with the NLL’ analytic calculation.

In this letter we perform the corresponding calculation for the LHCb experiment using two different methods. Our first method, which we will refer to as GFIP, is analogous to the GFIP calculation of Ref. Bain et al. (2016a). We start by generating events corresponding to hard production of $c$ quarks and gluons in $pp$ collisions at $\sqrt{s}=13$ GeV using MadGraph Alwall et al. (2014).222Contributions to $J/\psi$ production from quarks other than $c$ in the hard process are suppressed, either due to soft gluon emission or by $\alpha_{s}$ evaluated at a large energy scale. We therefore we neglect their contribution. In the LHCb data, all jets have pseudorapidity $2.5<\eta<4.0$, $R=0.5$, and the jets are required to have $p_{T}>20$ GeV. The hard partons generated by MadGraph satisfy the jet constraints of LHCb. PYTHIA is then used to shower the event down to a scale of $\sim 2m_{c}$. Jet algorithms are applied to the output of the PYTHIA shower and the $c$ quarks and gluons must be within jets of radius $R=0.5$ satisfying the criteria of the LHCb data described above. The resulting $c$ and gluon distributions are shown in Fig. 1. Note that the $c$ quark distribution is peaked near $z=1$ while the gluon $z$ distribution is much softer and peaked near $z=0$.

The $p_{T}$ and $y$ distributions for the $c$ quarks and gluons are then convolved manually with the NRQCD fragmentation functions evaluated at leading order (LO) in perturbation theory to obtain $p_{T}$ and $y$ distributions for $J/\psi$. For gluons we include ${}^{3}S_{1}^{[1]}$, ${}^{3}S_{1}^{[8]}$, ${}^{1}S_{0}^{[8]}$ and ${}^{3}P_{J}^{[8]}$ fragmentation functions, because the $v^{4}$ suppression of the color-octet LDMEs is compensated by powers of $\alpha_{s}$ for ${}^{1}S_{0}^{[8]}$ and ${}^{3}P_{J}^{[8]}$ and $\alpha_{s}^{2}$ for ${}^{3}S_{1}^{[8]}$. See Ref. Baumgart et al. (2014) for the explicit expressions for the LO NRQCD fragmentation functions. Color-singlet and color-octet fragmentation functions start at the same order in $\alpha_{s}$ for charm quarks so we include only color-singlet fragmentation for charm quarks. LHCb requires both muons have $2.0<\eta<4.5$, $p>5$ GeV, and $p_{T}>0.5$ GeV. The energy cut clearly suppresses contributions from partons with low $z$ and hence enhances the contribution from $c$ quark initiated jets. We implement the muon cuts by assuming the $J/\psi$ are unpolarized and therefore decays to $\mu^{+}\mu^{-}$ isotropically in its rest frame, and the LHCb cuts on the muons are applied to the muons after they are boosted back to the lab frame. From this a normalized distribution in $z(J/\psi)$ is constructed for each production mechanism. Each mechanism is characterized by an initial parton $i$ and quantum numbers $n$, and is multiplied by a weight

[TABLE]

The weight in Eq. (2) ensures that the total number of $J/\psi$ coming from each mechanism are in the proper ratio where $D_{i\to J/\psi}^{n}(z)$ are calculated at the scale $2m_{c}$. This is where the fitted LDMEs enter the calculation as $D_{i\to J/\psi}^{n}(z)\propto\langle{\cal O}^{J/\psi}(n)\rangle$. The LHCb data is normalized so that the sum of the heights of the bins adds to 1. Because of possible large corrections near $z\to 0$ and $z\to 1$, we only compare with LHCb data in the range $0.1<z<0.9$ and normalize our distributions to the sum of the data in these bins.

Our second method, which we refer to as the FJF method, employs FJFs combined with hard events generated by Madgraph at LO. In calculating the FJFs, logartihms of $m_{J/\psi}/E_{J}$ are resummed using leading order DGLAP equations to evolve the fragmentation functions from the scale $2m_{c}$ to the jet energy scale, $E_{J}$. Madgraph calculates the remaining terms in the factorization theorem to LO in perturbation theory. This does not include NLL’ resummation for the remaining terms in the factorization thereom, however the $z(J/\psi)$ dependence of the cross section is controlled primarily by the FJF. The energy distribution of hard partons is combined with the FJFs for anti-$k_{T}$ jets Waalewijn (2012) with $R=0.5$ to produce a $z(J/\psi)$ distribution for each of the five mechanisms. From the GFIP calculations, we know as a function of $z$ the fraction of $J/\psi$ that survive the muon cut and we apply this correction to our analytic calculations. The $z(J/\psi)$ distributions from each mechanism are weighted by the factors in Eq. (2) as before. The FJF is appropriate for $n$-jet cross sections like Eq. (1). Inclusive FJFs Kang et al. (2016, 2017a); Dai et al. (2016, 2017) differ by a contribution from out-of-jet radiation that is power suppressed for $R\sim O(1)$ Ellis et al. (2010).

Fig. 2 shows the predicted $z(J/\psi)$ distributions for the three choices of LDME’s in Table 1 using the GFIP (gray) and FJF (red) methods, which are in good agreement. Uncertainties are due to the LDMEs only. In the case of Ref. Bodwin et al. (2014), the errors in Table 1 are supplemented with an error correlation matrix Bodwin . In Ref. Chao et al. (2012) a fixed relationship between the ${}^{3}S_{1}^{[8]}$ and ${}^{3}P_{J}^{[8]}$ LDMEs is required to obtain unpolarized $J/\psi$. This constraint is taken into account when computing the uncertainty due to the LDMEs. These constraints significantly reduce the uncertainty in the predictions relative to naively adding uncertainties in Table 1 in quadrature. Other sources of uncertainty such as scale variation have not been included. Estimating theory uncertainties reliably in the absence of a complete factorization theorem is difficult. For example, using the FJF method, the $\mu$ dependence of the FJF should be cancelled by $\mu$ dependence in hard and soft functions that have not been computed. Note that since the normalization of theoretical curves is fixed to the LHCb data, any scale variation that affects normalization but not the shapes of the $z(J/\psi)$ distribution will not contribute to the uncertainty. Especially at low values of $z$, the underlying event and double parton scattering give additional theoretical uncertainties. However, it is not clear how estimate these uncertainties.

All three choices of LDMEs give better agreement to the LHCb data than default PYTHIA shown in Ref. Aaij et al. (2017). This gives support to the picture of quarkonium production in Ref. Bain et al. (2016a) and this letter. The LDMEs from global fits Butenschoen and Kniehl (2011, 2013) give worse agreement than the fits from Refs. Bodwin et al. (2014); Chao et al. (2012). The LHCb data is a decreasing function of $z(J/\psi)$ as $z(J/\psi)\to 1$. This is a property of the ${}^{3}S_{1}^{[1]}$ and ${}^{1}S_{0}^{[8]}$ FJFs, but not the ${}^{3}S_{1}^{[8]}$ and ${}^{3}P_{J}^{[8]}$ FJFs, which actually diverge as $z\to 1$. In order to obtain negligible polarization at high $p_{T}$, the ${}^{3}S_{1}^{[8]}$ and ${}^{3}P_{J}^{[8]}$ LDMEs of Refs. Bodwin et al. (2014); Chao et al. (2012) have relative signs such that they roughly cancel, so the ${}^{1}S_{0}^{[8]}$ dominates production and $J/\psi$ are unpolarized. The same cancellation here allows the $z(J/\psi)$ distribution go to zero as $z(J/\psi)\to 1$. Such a cancellation does not occur for the LDMEs from the global fits so the $z(J/\psi)$ distribution starts to turn up at large $z(J/\psi)$.

To summarize, we have analyzed the recent LHCb data on $J/\psi$ production within jets. We used a combination of Madgraph, PYTHIA, and LO NRQCD fragmentation functions first introduced in Ref. Bain et al. (2016a) as well as an approach based on Monte Carlo evaluation of the hard process combined with $J/\psi$ FJFs evaluated at the jet energy scale. Both methods yield $z(J/\psi)$ distributions that agree much better with data than default PYTHIA simulations. The $z(J/\psi)$ distributions are very well described by LDMEs from fits to large $p_{T}$ data, and less well described by LDMEs from global fits. It would be interesting to perform a combined fit to the LHCb data and the large $p_{T}$ data used in Refs. Bodwin et al. (2014); Chao et al. (2012) to see if consistent LDMEs with smaller errors can be obtained. Experimental measurement of jets at central rapidity and the polarization of $J/\psi$ as a function of $z(J/\psi)$ Kang et al. (2017b) would also be of interest. Finally it would be especially interesting to find ways of discriminating charm and gluon initiated jets Ilten et al. (2017), as a sample containing only gluon initiated jets will have greater sensitivity to color-octet LDMEs.