Probing Quarkonium Production Mechanisms with Jet Substructure

We use fragmenting jet functions (FJFs) in the context of quarkonia to study the production channels predicted by NRQCD (3S_1^(1), 3S_1^(8), 1S_0^(8), 3P_J^(8)). We choose a set of FJFs that give the probability to find a quarkonium with a given momentum fraction inside a cone-algorithm jet with fixed cone size and energy. This observable gives several lever arms that allow one to distinguish different production channels. In particular, we show that at fixed momentum fraction the individual production mechanisms have distinct behaviors as a function of the the jet energy. As a consequence of this fact, we arrive at the robust prediction that if the depolarizing 1S_0^(8) matrix element dominates, then the gluon FJF will diminish with increasing energy for fixed momentum fraction, z, and z>0.5.

(Dated: July 3, 2014) We use fragmenting jet functions (FJFs) in the context of quarkonia to study the production channels predicted by NRQCD ( J ). We choose a set of FJFs that give the probability to find a quarkonium with a given momentum fraction inside a cone-algorithm jet with fixed cone size and energy. This observable gives several lever arms that allow one to distinguish different production channels. In particular, we show that at fixed momentum fraction the individual production mechanisms have distinct behaviors as a function of the the jet energy. As a consequence of this fact, we arrive at the robust prediction that if the depolarizing 1 S (8) 0 matrix element dominates, then the gluon FJF will diminish with increasing energy for fixed momentum fraction, z, and z > 0.5. a Electronic address: baumgart@cmu.edu b Electronic address: akl2@pitt.edu c Electronic address: mehen@phy.duke.edu d Electronic address: izr@andrew.cmu.edu

I. INTRODUCTION
Nonrelativistic QCD (NRQCD) is an effective field theory [1] for quarkonium that reproduces full QCD as an expansion in the relative velocity, v, of the heavy quark and antiquark. This theory has been used to study both the decay and production of these bound states [2]. Its predictive power is predicated on our knowledge of a set of non-perturbative matrix elements that must be extracted from the data. In the case of J/ψ or Υ production there are four such matrix elements that must be fit at leading order, and thus predictions have mainly been limited to shapes of spectra. NRQCD predictions at NLO in the coupling have been compared to the world data on J/ψ production in Refs. [3,4]. The χ 2 /d.o.f. of 4.42 1 found in Ref. [3] is higher than one would hope for, but not unexpected given large theoretical uncertainties.
Thus, it is perhaps fair to say that we cannot yet claim that NRQCD is correctly describing quarkonium production with unqualified success. In particular, one prediction [5] of the theory is that, at asymptotically large transverse momentum, the 3 S 1 state (J/ψ or Υ) should be purely transverse at leading order. At present the data in both the charm and bottom sector do not see this trend [6] though the error bars are large, especially in the bottom system. Furthermore the various experiments are not in agreement.
It is important to appreciate that concluding that NRQCD is "wrong", in any sense, is equivalent to saying that QCD does not properly describe these states. If NRQCD predictions for large p T production are not agreeing with the data, and we assume that the data is correct, then the only logical alternatives are: (1) the velocity and/or α s expansions are not converging, (2) the fragmentation approximation, along with its expansion in m Q /p T , is wrong, either due to the failure of factorization or the presence of anomalously large power corrections. Let us consider each of these possibilities in turn. The perturbative corrections in α s (2m Q ) to the fragmentation function were found to be small [7]. The possibility that the velocity power counting could not apply to the charmed system [8,9] is certainly a viable option, though the velocity expansion seems to work relatively well for the decay processes [10]. Moreover, one would expect for the bottom system that the velocity expansion should converge nicely. It is possible that factorization is breaking down in the production processes, as all such proofs, at least within the confines of SCET, are lacking a treatment of the factorization breaking "Glauber mode". Nonetheless, given the success of semi-inclusive predictions in light hadronic systems, it would be surprising to see a failure in the case of quarkonium.
A more conservative guess would be that there is nothing wrong with the theory, but perhaps the values of extracted matrix elements are sufficiently inaccurate as to change the nature of the polarization prediction. For example, the magnetic spin flip operator could be anomalously large. In any case, to get a better handle on the situation we must improve our quantitative understanding of the various production channels associated with the aforementioned matrix elements. The purpose of this paper is to introduce a new tool that will allow for a new extraction of these matrix elements by studying the characteristics of jets within which the quarkonium reside.

II. THE FRAGMENTING JET FUNCTION (FJF)
Power counting dictates that at asymptotic values for p ⊥ m Q , quarkonia should be produced by single parton fragmentation. 2 Since the parton initiating the fragmentation is a colored object, the quarkonium will be produced in association with light hadrons. In this paper we will consider a J/ψ 3 produced within a jet of energy E and cone size R, in which the J/ψ carries a fraction of the jet energy, z. In this situation, a generic cross section is determined by the convolution of a hard and soft function (and possibly other jet functions, if there are other jets detected in the final state) multiplied by a quantity known as the fragmenting jet function (FJF), first introduced in Ref. [14] and further studied in Refs. [15][16][17][18][19]. These papers focused on FJFs for light hadrons such as pions. FJFs for particles with a single heavy quark are studied in Ref. [20]. We show that the FJFs for gluon and charm quark jets containing a J/ψ can be calculated in terms of set of NRQCD long-distance matrix elements (LDME). In our calculations the relevant LDMEs are: 0 ) . The spectroscopic notation indicates the quantum numbers of the heavy quarks prior to hadronization. We show that the contribution to the FJF from each of these mechanisms depends differently on z and E and can thus be used to extract the LDME. Our results could easily be extended to jets containing other quarkonia states.
Since there are many observables associated with jets (angularities [21], broadening [22], jet shape [23], Nsubjettiness [24], etc.), one can generate a very large number of new tests of the NRQCD factorization formalism by applying jet physics techniques to the study of quarkonia produced within jets. Furthermore, studying high p ⊥ quarkonia produced within jets avoids some of the potential theoretical pitfalls that could plague tests of the NRQCD factorization formalism at small p ⊥ . At the highest p ⊥ available, we expect factorization to hold up to corrections which scale as m Q /p ⊥ , and furthermore the α s expansion should be well behaved.

A. Operator Definitions
We first briefly review the properties of the FJF [14][15][16][17][18][19]. While we will not explicitly work in the effective field theory SCET [25], we will always have in mind the idea that one can construct factorization theorems within this formalism. The factorization theorem for the production cross section for a jet with energy E, cone size R, and a J/ψ with energy fraction z in a pp collision is schematically of the form where H is the hard function, f a/p and f b/p are parton distributions functions, J i are jet functions for each jet not containing the J/ψ, G ψ is the FJF for the jet containing the J/ψ, and S is the soft function. Generically there are two types of jets, unmeasured and measured, in the terminology of Ref. [26]. Unmeasured jet functions describe jets in which only the large light-cone momentum (measured along the jet axis) is known. In measured jets, some aspect of the jet's substructure has also been measured. For unmeasured jets, soft gluon radiation does not affect the total momentum of the jet (up to power corrections) and therefore these jet functions enter the cross section multiplicatively. For measured jets, the jet substructure may be sensitive to the soft radiation, therefore it must be convolved with the soft function. For G ψ (E, R, z, µ), R, E and z are not affected by soft radiation (up to power corrections) so it enters the cross section multiplicatively and all of the z dependence is contained in G ψ (E, R, z, µ), which enables us to ignore all the other factors in Eq. (1) and focus on G ψ (E, R, z, µ).
A generic fragmenting jet function may be defined as a product of operators of the form where O int is some interpolating field for the parton of interest. O meas is a measurement operator (a set of delta functions) that fixes the measured jet characteristics, such as E, R and z. The operators are manifestly gauge invariant. In SCET these operators would involve only fields with the same large momentum (and possibly soft fields) and compose a piece of the factorization theorem (not shown) that is generated at the highest scale Q, which is usually taken to be on the order of the jet energies E. 4 G ψ i contains two relevant scales: the invariant mass or energy of the jet and the hadron mass. Thus one can perform a further factorization to separate out these two scales where the long distance physics is captured by a fragmentation function, and the short distance physics can be calculated perturbatively. The resulting form of this second step of factorization can be written as [14] where we have now specialized to the case of interest where the jet energy E is measured for cone size R. Loosely speaking this function gives the probability of finding a quarkonium whose large momentum fraction, relative to the jet within which it is found, is z. The operator definition of the quark fragmentation function is [27] given by where the operator q includes an anti-path ordered Wilson line that renders the matrix element gauge invariant. A similar matrix element can be written down for the gluon fragmentation function. What distinguishes the quarkonium fragmentation function from other cases is that it contains a further subset of scales: the quark mass, the Bohr radius, and the binding energy that scale as 1, v, and v 2 respectively in units of the quark mass. Furthermore, taking the quark mass scale to be perturbative implies that the constituents are produced at a point, and that the momentum fraction carried by the quarkonium is set perturbatively. This is so even if the pair is produced in an octet state, since the shedding of color occurs via soft multipole emission whose effect on the kinematics is suppressed by an amount of order v 2 , except near the end point z = 1 where these non-perturbative corrections are enhanced and can be accounted for by the inclusion of a non-perturbative shape function [28]. In general we will present our results away from the end point to avoid the need for such a function. Thus, the fragmentation functions for quarkonium are calculable up to a set of LDMEs. The matching coefficients J ij (E, R, z, µ) can be calculated in perturbation theory. Large logarithms in the J ij (E, R, z, µ) are minimized at the scale 2E tan(R/2)(1 − z). Note that the matching coefficients J ij (E, R, z, µ) are independent of the choice of hadronic final states, and thus we may utilize the results in Ref. [18] for the FJF for light hadrons for the case at hand. B. Expressions for the J/ψ FJF We will focus gluon and charm quark fragmentation to J/ψ. For gluon fragmentation to J/ψ through cc pairs, we consider the 3 S gluon fragmentation function is leading order in the v expansion, as the color-octet contributions are suppressed by v 4 . However the gluon color-singlet contribution is suppressed relative to 3 S contributions because both color-singlet and color-octet mechanisms start at the same order in α s . The ratio of gluon to charm production cross sections at the LHC is approximately 50, but the ratio of charm quark to gluon fragmentation functions, partially compensates for this suppression. Fragmentation from light quarks is suppressed by one power of α s relative to the 3 S gluon fragmentation contribution and shares the octet velocity suppression. The J ij (E, R, z, µ) and the relevant fragmentation functions are collected in the Appendix.
The convolution in Eq. (3) can be explicitly evaluated using the formula for J gg (E, R, z, µ) and J gq (E, R, z, µ) in the Appendix to obtain +θ In this expression, we have defined This expression shows that the logarithms in G ψ g (E, R, z, µ) are minimized at the scale µ = 2E tan(R/2)(1 − z), as first pointed out in Ref. [18]. The logarithms of 1 − z are easily resummed using the jet anomalous dimension [18], however, we will not do this resummation in this paper as we consider 1 − z ∼ O(1). We instead set the scale in J gg (E, R, z, µ) to be µ J = 2E tan(R/2), and evolve the fragmentation function from the scale 2m c to the scale µ J . This is done by taking moments of the fragmentation functions, evolving each moment according to its anomalous dimension as obtained from the Altarelli-Parisi equations, and then performing an inverse-Mellin transform.
The G ψ g(q) (E, R, z, µ) is present because of mixing with the quark fragmentation function. In principle there should be a sum over all quark flavors. However, the light quark fragmentation function contributes only via fragmentation through 3 S gluon fragmentation so it will be neglected. Charm quarks and antiquarks can fragment via 3 S (1) 1 cc pairs at O(α 2 s ), which is lower order than the corresponding gluon fragmentation function. Therefore this mixing must be included.
The quark FJF is given by: For this contribution, as previously mentioned we will only consider the q = c contribution fragmenting via 3 S (1) 1 cc pairs. The mixing contribution of gluon fragmentation into this FJF must also be included. To evaluate G ψ i (E, R, z, µ J ) we will use Eqs. (5-8) with our numerically evaluated D i→ψ (z, µ J ). We see that up to O(α s ) corrections which shows that the z distribution of a J/ψ within a jet with energy E and cone size R is approximately equal to the fragmentation function evaluated at the jet scale µ J = 2E tan(R/2). Since the fragmentation functions for 3 S are very different, this observable has the power to discriminate between all four gluon-production mechanisms. This can seen from a cursory inspection of the expressions for the fragmentation functions given in the Appendix and shown in Fig. 1. Though the dramatic differences in these functions are considerably softened by Altarelli-Parisi evolution, we will see that each contribution to G ψ g (E, R, z, µ) has a different E dependence that varies for fixed z (cf. Fig. 3). This makes it clear that measurement of G ψ g (E, R, z, µ) for different momentum fractions has potential to allow independent extraction of all four LDME. In our calculations E and R will always enter in the combination E tan(R/2) and we will choose R = 0.4 .
In Fig. 2 we plot the 3 S         It is also interesting to study the energy dependence of the fragmentation functions. In Fig. 3 we plot the four gluon FJFs as a function of energy E for three different values of z using the same color-coding as above. The LDME of Refs. [3,4] have again been used to set the normalization of the curves. In order to the make shapes of the curves more easily viewable, we have divided the 3 P  is decreasing only for 0.8. Extractions of the E dependence of the FJF for different values of z should allow one to disentangle the various contributions to quarkonium production. In particular, note that if the lack of polarization is due to an anomalously large 1 S (8) 0 , then we should see a decrease in the gluon FJF as a function of the jet energy for fixed z, with z > 0.5.
The moments of the FJF, z N ≡ 1 0 dz z N −1 G ψ g (E, R, z, µ), can be calculated analytically using the formulae in the Appendix. Note that this integral diverges if N = 1 because the N = 1 moments of both the Altarelli-Parisi splitting function and the matching coefficients J gg (E, R, z, µ) have poles at N = 1. This could be cured by resummation of log z, as implemented for the D g→ψ (z, µ) fragmentation function in Ref. [29], but this is beyond the scope of this paper. The LDME cancel in the ratios of moments, and we plot ratios of successive moments, z N +1 / z N , for N = 2, 3, and 4 in Fig. 4. In all columns we have plotted the moment ratios of the 3 S (1) 1 FJF (black). We also plot moment ratios for the 3 S  moments have power to discriminate between various production mechanisms, in particular, we find but once scale uncertainties are included it is hard to distinguish these two moment ratios. The energy dependence of the moments of the color-octet FJFs is given by When we set µ ≈ 2E tan(R/2), the energy dependence is entirely contained in the first two factors on the r.h.s. of Eq. (12), which are the same for all three color-octet FJFs. The color-singlet and charm quark fragmentation functions are more complicated due to the mixing of these fragmentation function in the evolution from the scale 2m c to µ J . Making log-log plots of z N we find that that z N ∝ (log E) F (N ) where F (N ) can be extracted from Eq. (12).

III. COMPARISON OF VARIOUS LDME EXTRACTIONS
In the final part of this paper, we will discuss what recent extractions of the LDME predict for the gluon FJF. In addition to the extractions in Refs. [3,4], we will consider values of the LDME extracted in two recent papers [30,31] that attempt to solve the polarization puzzle by focusing exclusively on high p ⊥ production of charmonia at collider experiments. The study in Ref. [30] uses a NLO NRQCD calculation to fit the color-octet LDME to inclusive J/ψ production at high p ⊥ and finds values of the LDME that can produce negligible polarization in agreement with the data. However, these values of LDME are inconsistent with the results of fitting the world data in Refs. [3,4]. In particular, O J/ψ ( 1 S (8) 0 ) is larger by a factor of two and O J/ψ ( 3 P (8) 0 ) has the opposite sign as the fit in Refs. [3,4]. These two effects combine to produce significant depolarization of the J/ψ. In Ref. [31], the calculations are performed in the leading-power fragmentation approximation and logarithms of p ⊥ /m c are resummed by using Altarelli-Parisi equations for the fragmentation functions. The fitted LDME are similar to those found in Ref. [30] in the sense that O J/ψ ( 1 S 1 . It should be noted that the quoted errors in the extracted LDME in Refs. [3,4] are considerably smaller than those in Refs. [30,31]. However, the presence of nontrivial correlations between the uncertainties in [31] allows us to make a much sharper prediction for the gluon FJF than is naively suggested by the large individual error bars [32]. In all of these extractions, there is a hierarchy between matrix elements that are supposed to have the same velocity scaling. However, it is generated by anomalously small matrix elements not anomalously large ones. In Fig. 5, we compare the predictions for the gluon FJF at E = 50 GeV and E = 200 GeV using the results from the fits to the LDME in Refs. [3,4,30,31]. The gluon FJF is the sum over all contributions, color-singlet as well as color-octet. The color-singlet matrix element is chosen to be 1.32 GeV 3 in Refs. [3,4,31] and 1.16 GeV 3 in Ref. [30]. We use the LDME extracted in the original fit and the error bands are the result of adding in quadrature the uncertainties for the LDME quoted in Refs. [3,4,30]. We supplement the uncertainty given in Ref. [31] with the full correlation matrix provided by one of the authors [32]. No other theoretical uncertainty is included. The gray The gluon FJF at fixed momentum fraction for the LDME extracted in Refs. [3,4] (gray), Ref. [30] (blue), and Ref. [31] (red).
band with black borders is the prediction using the LDME extracted in Refs. [3,4], the red band uses the matrix elements extracted in Ref. [31] and the blue band uses the matrix elements extracted in Ref. [30]. Fig. 6 shows the energy dependence at fixed momentum fraction for the different determinations. We see that for z >0.5, the question of which set of LDMEs is preferred, those determined for the world average [3,4] or those that alleviate the polarization puzzle [30,31], will be resolved by testing whether the gluon FJF is increasing or decreasing with energy. Furthermore, measurement of the gluon FJF has the power to distinguish between all three fits.

IV. CONCLUSIONS
We have demonstrated that by studying the characteristics of jets arising from quarkonium production, we can disentangle the various production channels. There are a multitude of ways of analyzing such events. Here we have chosen to measure the energy and cone angle of the jet, but one could consider other observables such as the invariant mass. Within our choice of variables (E, R) we found that a particularly discriminating tool is the measurement of the energy dependence at fixed momentum fraction as shown in Figs. 3 and 6. A robust prediction of our analysis is that for z > 0.5 the gluon FJF at fixed z should decrease as function of energy if the lack of transverse polarization in the data is due to the dominance of the 1 S (8) 0 LDME over the other color octet matrix elements for high-p ⊥ production. Further information can be gathered by calculating the normalized cross section, in which case one could constrain the sum of the matrix elements.
The integrals over r and y must be done numerically. The 1 S (8) 0 gluon fragmentation function is given by [35][36][37] and the 3 P gluon fragmentation function is given by Here we have summed over J = 0, 1, 2 and used O ψ ( 3 P charm quark fragmentation function is [34], The moments of the color-octet gluon fragmentation functions can be computed analytically. Defining we haveD The fragmentation function is evolved using the standard DGLAP evolution, µ ∂ ∂µ D i (z, µ) = α s (µ) π j 1 z dy y P i→j (z/y, µ)D j (y, µ) .
These equations are solved analytically in moment space, and then the fragmentation functions at the scale µ J are obtained by numerically evaluating the inverse Mellin transform. In our calculations, q = c and the mixing between the gluon and c quark fragmentation function is only relevant for the 3 S (1) 1 channel.
It is useful to have analytic expressions for the moments of the matching coefficients; these are given by: J gg (E, R, N, µ) = 1 0 dz z N −1 J gg (E, R, z, µ) 2(2π) 3 = 1 + α s C A π L 2 + P N gg L + H 2 N −1 − 5π 2 24 + H N −1,2 J qq (E, R, N, µ) = 1 0 dz z N −1 J qq (E, R, z, µ) 2(2π) 3 = 1 + α s C F π L 2 + P N qq L + where H N is the harmonic number, H N,2 is the generalized harmonic number of order 2, and P N ij , F N , and G N are given by .
Note that so in the large N limitJ gg (E, R, N, µ) = 1 + J qq (E, R, N, µ) = 1 + α s C F π L 2 N − where L N = ln 2E tan(R/2) N e γ E µ . (A32) We see that the logarithms inJ gg (E, R, N, µ) andJ qq (E, R, N, µ) are minimized at the scale 2E tan(R/2)/N e γ E . This is consistent with the expressions for J gg (E, R, z, µ) and J qq (E, R, z, µ) These logarithms are easily resummed using the jet anomalous dimension, however, we will not do this resummation in this paper as we compute moments with N of order unity. Moments with N = 1 are divergent because of the poles in P N gg , P N qg ,J gg , andJ gq .