Inclusive $J/\psi$ and $\eta_c$ production in $\Upsilon$ decay at $\mathcal{O}(\alpha_s^5)$ in nonrelativistic QCD factorization

We study $J/\psi$ and $\eta_c$ inclusive production in $\Upsilon$ decay within the framework of nonrelativistic-QCD (NRQCD) factorization. In the latter case, for which no experimental data exist so far, we also include the $h_c$ feed-down contribution. We calculate the short distance coefficients completely through $\mathcal{O}(\alpha_s^5)$. The NRQCD predictions for the branching fraction $\mathcal{B}(\Upsilon\to J/\psi+X)$ via direct production, evaluated with different sets of long-distance matrix elements (LDMEs), all agree with the experimental data in a reasonable range of renormalization scale. Using $\eta_c$ and $h_c$ LDMEs obtained from $J/\psi$ and $\chi_c$ ones via heavy-quark spin symmetry, we find that the bulk of $\mathcal{B}(\Upsilon\to\eta_c+X)$ via prompt production arises from the $c\bar{c}({}^3\!S_1^{[8]})$ Fock state. The experimental study of this decay process would, therefore, provide a particularly clean probe of the color octet mechanism of heavy-quarkonium production.


I. INTRODUCTION
Heavy-quarkonium production serves as an ideal laboratory to study the interplay of perturbative and nonperturbative phenomena of QCD thanks to the hierarchy of energy scales m Q v 2 Q ≪ m Q v Q ≪ m Q characterizing kinetic energy, momentum, and mass, where m Q is the mass of the heavy quark Q and v Q is its relative velocity in the rest frame of the heavy quarkonium. The effective quantum field theory of nonrelativistic QCD (NRQCD) [1] endowed with the factorization conjecture of Ref. [2] is the standard theoretical approach to study quarkonium production and decay. This conjecture states that the theoretical predictions can be separated into processdependent short-distance coefficients (SDCs) calculated perturbatively as expansions in the strong-coupling constant α s and supposedly universal long-distance matrix elements (LDMEs), scaling with definite powers of v Q [3]. In this way, the theoretical calculations are organized as double expansion in α s and v Q .
During the past quarter of a century, the NRQCD factorization approach has been very successful in describing both heavy-quarkonium production and decay; see Refs. [4][5][6] for reviews. However, there are still open problems in charmonium production, in particular for the J/ψ meson. Prompt J/ψ production has been studied in various scenarios both experimentally and theoretically. Specifically, the SDCs are known at next-to-leading order (NLO) in α s for the yield [7,8] and polarization [9] in e + e − annihilation, the yield in two-photon collisions [10][11][12], the yield [13] and polarization [14] in pho- * Electronic address: zhiguo.he@desy.de † Electronic address: kniehl@desy.de ‡ Electronic address: xiangpeng.wang@anl.gov toproduction, the yield [15,16] and polarization [17][18][19][20] in hadroproduction, etc. Different sets of LDMEs were obtained by fitting experimental data adopting different strategies. Unfortunately, none of them can explain all the experimental measurements, which challenges the universality of the NRQCD LDMEs. Moreover, it has been found [21] that all the LDME sets determined from J/ψ production data, upon conversion to η c LDMEs via heavy-quark spin symmetry, result in NLO predictions that overshoot the η c hadroproduction data, taken by the LHCb Collaboration at the CERN Large Hadron Collider (LHC) [22]. To shed more light on this notorious J/ψ puzzle, it is useful to consider yet further production modes. This provides an excellent motivation to study J/ψ and η c inclusive production in Υ decays, which is the goal of this paper. This complements our recent study of χ cJ inclusive production in Υ decays [23]. To facilitate the comparison with experimental data, we will focus here on direct J/ψ production, while we will consider prompt η c production, by also including the feed-down from h c mesons.
On the experimental side, the Υ → J/ψ+X decay process was first observed by the CLEO Collaboration at the Cornell Electron Storage Ring CESR about three decades ago [24]. Several independent measurements [25][26][27] followed. The latest one was carried out by the Belle Collaboration at the KEK B factory KEKB [28]. The present world average of the branching fraction for prompt production is [29] B prompt (Υ → J/ψ + X) = (5.4 ± 0.4) × 10 −4 . (1) As for the Υ → η c + X decay, there is no experimental data yet. With the large amount of data to be accumulated by the Belle II Collaboration at the SuperKEKB accelerator, there might be a chance to study this process, and NRQCD predictions for its branching fraction will be needed.
On the theoretical side, the Υ → J/ψ + X process was proposed as a rich gluon environment to study the color octet (CO) mechanism of NRQCD factorization in Refs. [30,31], where some of the CO channels, including bb( 3 S s ), were taken into account. Later, the color singlet (CS) processes, including bb( 3 S s ), and some interesting QED contribution were also calculated [32,33]. It was found that the CS contribution itself is about 3.8 times smaller than the experimental data, which suggests the potential need for a large CO contribution.
In the previous studies [30][31][32][33], only the leading-order (LO) O(α 4 s ) contribution and some part of the O(α 5 s ) channels were considered and pre-LHC LDME sets were used. In this work, we will obtain the complete SDC results through O(α 5 s ), both for J/ψ and η c production, and perform a numerical analysis with up-to-date LDMEs sets. On top of direct production, prompt production also includes the feed-down from heavier charmonia, namely from χ cJ and ψ ′ mesons in the J/ψ case and from h c mesons in the η c case. The direct J/ψ results obtained here readily carry over to ψ ′ feed-down, while the χ cJ results may be found in Ref. [23]. The h c feed-down results will also be derived here.
This paper is organized as follows. In Sec. II, we describe our analytical calculations. In Sec. III, we present and interpret our numerical results. Section IV contains our conclusions. In the Appendix, we list the contributing SDCs at O(α 3 s ).

II. ANALYTICAL CALCULATIONS
According to NRQCD factorization [2], the partial decay width of Υ → H + X with H = J/ψ, η c , h c is given by whereΓ mn =Γ(bb(m) → cc(n) + X) is the SDC and Υ|O(m)|Υ and O H (n) are the LDMEs of Υ decay and H production, respectively.
For Υ, the CO contribution is so small [30] that we only include the CS case of m = 3 S 1 . In fact, by the velocity scaling rules [3], the leading CO LDMEs of Υ decay are parametrically suppressed by a factor of v 4 b ≈ 1% relative to the CS one, Υ|O( 3 S 1 )|Υ . In Ref. [34], the suppression factors for the S wave CO LDMEs of Υ decay were quantitatively determined in the framework of lattice NRQCD using two different schemes of implementing heavy-quark Green's functions, the hybrid and nrqcd schemes. The largest values were obtained in the hybrid scheme, namely, Υ|O( 1 S  Table I in Ref. [34]. Similar analyses for the P wave LDME Υ|O( 3 P [8] 0 )|Υ are not yet available, but its suppression factor is expected to be in the same ballpark, of O(10 −3 ) or below. On the other hand, SDCs involving CO bb Fock states m already start at O(α 3 s ). However, using α s (m b ) ≈ 0.26, we have α 2 s ≈ 0.068 ≫ 10 −3 , so that the contributions to Eq. (2) from the leading CO bb Fock states, m = 1 S [8] 0 , 3 S 1 , which we still include here. For future use, we list the O(α 3 s ) SDCsΓ mn for all the decay processes considered here and in Ref. [23], On the other hand, for H, we include the CO contributions of LO in v 2 c besides the CS contributions, i.e. we have n = for η c , and n = 1 P 0 for h c . In our computation, we generate all the Feynman amplitudes using the FeynArts package [35]. The Dirac and color operations are performed with the help of the Feyn-Calc [36] and FORM [37] packages. We use the Mathematica package $Apart [38] to decompose linearly dependent propagators in the loop integrals to irreducible ones, which we further reduce to master scalar integrals using the FIRE package [39]. We then evaluate the master integrals numerically using the C++ package QCD-LOOP [40]. Finally, we perform the phase space integrations numerically using the CUBA [41] library.
The O(α 4 s ) LO contribution only includes the bb( 3 S 1 ) + gg subprocess. Its NLO QCD corrections, of O(α 5 s ), were calculated in our previous work [23], where details may be found. Here, we only present numerical results. At O(α 5 s ), also a number of new partonic production channels open up. We derive here the SDCs of all of them. In the cases where results already exist in the literature, we compare and find agreement. For the sake of a systematical discussion, we divide the contributing partonic subprocesses into three groups according to the parton content of the system X in the final state: (a) X = g, (b) X = ggg, and (c) X = ccg. Representative Feynman diagrams for each group are shown in Fig. 1.
In group (a), the only possible cc Fock states are n = 1 S J , and there are 6 Feynman diagrams, a typical one of which is shown in panel (a) of Fig. 1. The SDCs of these partonic subprocesses were first calculated in Ref. [31]. We evaluate the one-loop amplitude of bb( 3 S 0 ) + g applying the helicity projector approach [42] and the one of bb( 3 S J ) + g as the virtual corrections in Ref. [23]. We find agreement with Ref. [31].
In group (b), the possible cc Fock states are n = infrared finite and straightforwardly calculated. The infrared divergence in the SDC of bb( 3 S J )+ggg can be extracted by means of the phase space slicing method [43], as explained in detail for the CS case of J in Ref. [23]. In NRQCD factorization, this infrared divergence is absorbed into the NLO corrections to the where 1 ) ren is the renormalized LDME, µ and µ Λ are the QCD and NRQCD factorization scales, respectively, ǫ = 2 − d/2 in d space-time dimensions, and C F = 4/3 and C A = 3 are color factors. We verified numerically that the finite piece left over after removing the infrared divergence does not depend on the slicing parameter δ s , which is illustrated in Fig. 2.
In group (c), all the relevant cc Fock states n = J appear. There are generally 6 Feynman diagrams similar to the one shown in panel (c1) of Fig. 1, with the exception of n = 3 S 1 ) + ccg was first considered in Ref. [45], and the correct result was first presented in Ref. [32]. All the other partonic subprocesses are studied here for the first time. Their SDCs are all infrared finite and can be calculated straightforwardly.
To be able to investigate the relative weight of each partonic subprocess, we decomposeΓ 3 S [1] 1 ,n for each n according to the possible final-state systems X,  where X = gg, ggg, ccg for n = 3 S 1 . Furthermore, we factor out the strong-coupling constant α s (µ) and render the coefficients dimensionless by writinĝ where k = 5 for n = 3 S Given the large powers of α s (µ) appearing in Eq. (5), we are faced with considerable µ dependencies, so that scale optimization appears appropriate. As in Ref. [23], we adopt the principle of fastest apparent convergence (FAC) [46] to determine the default value of µ, by requiring that the LO and NLO evaluations ofΓ gg 3 S [8] 1 as described above coincide. This yields µ FAC = 6.2 GeV. It would then be natural to explore the µ dependence in the range µ FAC /2 < µ < 2µ FAC . However, the NLO result for B(Υ → χ c0 + X) is known to be negative for µ < 3.7 GeV [23]. We thus consider the reduced µ range 3.7 GeV < µ < 2µ FAC in the following.
When the bb pair is in a CS Fock state, it couples to at least three gluons in the decays Υ → H + X with H = J/ψ, η c , h c , as may be seen in Fig. 1. Therefore, the large theoretical uncertainty due to the freedom in the choices of α s and m b can be greatly reduced by normalizing Γ(Υ → H + X) with respect to Γ(Υ → ggg) [23,32,33]. We thus write the branching fractions of interest as where B(Υ → ggg) = 81.7% [29]. Through O(α 4 s ), we have [47] where β 0 = 11 − 2n f /3 with n f = 4. As an additional benefit, the theoretical predictions do not depend on the CS LDME Υ|O( 3 S 1 )|Υ anymore.
We first discuss Υ → J/ψ + X via direct production. From Table I J ) + ggg channels, the latter one being studied here for the first time. Several J/ψ LDME sets have been determined, both at LO and NLO, by fitting to different sets of J/ψ yield and polarization data. There are significant differences between them, which reflects the potential challenge to NRQCD factorization mentioned above. We select here the ones most frequently used in the recent literature [12,19,[48][49][50] and summarize their values in Table II. Besides the contributions through O(α 5 s ) discussed above, we will also include some other non-negligible contributions previously studied, namely the pure CS contributions at O(α 2 s α 2 ) [32] and O(α 6 s ) [33]. From Ref. [29], we extract the result for direct J/ψ production, by subtracting the contributions due to the feed-down from the χ cJ and ψ ′ mesons. In Fig. 3, we confront it with our NRQCD predictions evaluated with the J/ψ LDME sets listed in Table II as functions of µ. We observe from Fig. 3 that, although the various predictions appreciably differ in normalization and line shape, each of them can describe the experimental data in a plausible region of µ. For each of the considered J/ψ LDME sets, we determine the value of µ, µ fit , where the respective theoretical prediction coincides with the central value of the experimental result and list it in Table III. We observe from Table III that these µ fit values all fall into the preferred range 3.7 GeV < µ < 2µ FAC and thus conclude that NRQCD factorization is compatible with experiment here. However, we caution the reader that the J/ψ LDME sets of Refs. [12,19,49,50] yield µ fit values that are rather close to the border at 3.7 GeV, below which the NLO prediction for B(Υ → χ c0 + X) is negative [23].
We now discuss Υ → η c + X via prompt production. The LDMEs of η c and h c production are related to those of J/ψ and χ c production via heavy-quark spin symmetry We now evaluate B(Υ → η c + X) via direct production using the SDCs in Table I together with the η c LDMEs converted from Table II via Eq. (9) and present the results as functions of µ in Fig. 4. The corresponding results at µ = µ FAC are listed in Table IV. For comparison, we also present there the results at µ = µ fit for the individual µ fit values in Table III. We observe from Fig. 4 and the entries in Table IV for the common scale choice µ = µ FAC that the results for the J/ψ LDME sets from Refs. [19,50] and the default one from Ref. [18] are very similar, quite in contrast to the findings in Fig. 3 and Table III. This may be traced to the approximate agreement of the respective values of O J/ψ ( 1 S [8] 0 ) in Table II by noticing that, similarly to the case of η c hadroproduction [21], the bulk of the NRQCD prediction is made up by the 3 S 1 channel of η c production is quenched leading to a particularly small NRQCD prediction for B direct (Υ → η c + X). In the case of η c inclusive hadroproduction, it was found at NLO that the CS channel alone exhausts the cross section measured by the LHCb Collaboration [22], while the complete NRQCD predictions overshoot it by a few times to an order of magnitude depending on the LDME set chosen [21]. It is, therefore, interesting to examine the importance of CO processes for Υ → η c + X. We find that the CS contribution to B direct (Υ → η c + X) ranges from 3.2 × 10 −5 to 0.5 × 10 −5 for 3.7 GeV < µ < 2µ FAC and takes the value 1.1×10 −5 for µ = µ FAC , which is about 1.6 times smaller than the smallest NRQCD prediction for µ = µ FAC in Table IV, for Set 2 of the J/ψ LDMEs from Ref. [18], and more than one or even two orders of magnitude smaller than NRQCD predictions for the other J/ψ LDMEs sets. In other words, Υ → η c +X is predicted to proceed even more dominantly through CO processes than η c inclusive hadroproduction. We thus conclude that the study of Υ → η c + X is expected to provide a crucial test of the CO mechanism of NRQCD factorization and to place valuable constraints on the CO LDMEs of J/ψ production, in particular on O J/ψ ( 1 S [8] 0 ) , assuming that HQSS is preserved. Finally, we caution the reader that, in the case of B direct (Υ → η c + X), the CO contributions from the initial state may be somewhat less suppressed than expected from the velocity scaling rules [3]. In fact, comparing the entry for the J/ψ LDMEs of Refs. [12,49] and µ = µ FAC (µ fit ) in Table IV with the respective entry for the hybrid scheme in Table V, we observe that the contribution due to the S wave CO bb Fock states reaches about 26% (16%) of our default result. D. Υ → hc + X Finally, we turn to Υ → η c + X via the feed-down from h c mesons, i.e. via Υ → h c + X followed by h c → η c + γ. Converting the LDMEs for χ c0 production quoted above to those of h c production via the appropriate HQSS relations in Eq. (9), we obtain O hc ( 1 P sented as a function of µ in Fig. 6. At µ = µ FAC , we have B(Υ → h c + X) = 1.6 × 10 −5 . Taking into account the branching fraction B(h c → η c + X) = 51% [29], the feed-down contribution is found to be B hcfeed down (Υ → η c + X) = 8.0 × 10 −6 at µ = µ FAC .

IV. CONCLUSIONS
In summary, we studied J/ψ and η c inclusive production via Υ decay in the NRQCD factorization approach, including also the feed-down from h c mesons in the latter case. We calculated the SDCs completely through O(α 5 s ) keeping the bb pair in the 3 S  As for Υ → J/ψ + X via direct production, we found that the NLO sets of J/ψ LDMEs recently determined from fits to experimental data of inclusive production at various types of colliders or just hadron colliders all lead to useful interpretations of the measured branching fraction, in the sense that the renormalization scale µ fit which yields agreement with the central value is in the ballpark of the optimal scale µ FAC = 6.2 GeV of the FAC principle.
In want of direct determinations, we relied on η c and h c LDMEs obtained from J/ψ and χ cJ ones via HQSS relations. As for Υ → η c + X via prompt production, we thus predicted the branching fraction to range between 2.6 × 10 −5 and 2.2 × 10 −3 at µ = µ FAC , with a strong dependence on the choice of LDMEs, in particular on the value of O ηc ( 3 S [8] 1 ) . On the contrary, the CS contribution to B prompt (Υ → η c + X) amounts to just about 1.1 × 10 −5 at µ = µ FAC . Our findings suggest that inclusive η c production in Υ decay will provide a particularly clean laboratory to investigate the CO mechanism in charmonium production and so to shed light on the J/ψ polarization puzzle at hadron colliders, and we propose to carry out such measurements in the future.

Acknowledgments
We would like to thank G. T. Bodwin  The U.S. Government retains for itself, and others acting on its behalf, a paid-up nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government.
Appendix: SDCsΓmn at O(α 3 s ) Here we list our analytic results for the SDCsΓ mn of the decay processes Υ → H+X with H = J/ψ, η c , χ cJ , h c at O(α 3 s ). We havê Γ(bb( 1 S , Γ(bb( 3 S . (A.1) We note that, incidentally,Γ(bb( 3 S  from Refs. [12,49] and the χ cJ ones from Ref. [23] and deriving from them the η c and h c ones via the HQSS relations in Eq. (9), we find that the results for the ratios Υ|O( 1 S  Table V. These results are µ independent thanks to the normalization in Eq. (6), which, for consistency, is evaluated here at LO. They still need to be supplemented by the contributions proportional to Υ|O( 3 P