Abstract

In this paper, we review the current status of the phenomenological study of quarkonium production in high-energy collisions. After a brief introduction of several important models and effective field theories for quarkonium production, we discuss the comparisons between theoretical predictions and experimental measurements.

1. Introduction

Since the discovery of the in 1974, heavy quarkonium has been on the focus of much experimental and theoretical attention. Heavy quarkonium is a bound state consisting of a heavy quark () and its antiquark (). Depending on the flavor of the quark pair, there are charmonium and bottomonium. The production of a heavy quarkonium involves three different momentum scales: the heavy quark mass ( and in on-shell scheme), which governs the perturbative creation of the heavy quark pair (); the heavy quark momentum in the quarkonium rest frame; and the typical heavy quark kinetic energy , which governs the nonperturbative hadronization of the to physical quarkonium. Here, is the typical heavy quark velocity in the quarkonium rest frame ( for charmonium and for bottomonium). Due to the nonrelativistic nature of the bound state, heavy quarkonium production at high-energy collisions is a very important process to test our understanding of QCD.

Experimentally, many quarkonium states are relatively simple to identify in different colliders because of their clean experimental signatures and reasonably high yields. By virtue of these advantages, the heavy quarkonium is considered a promising tool to study the inner parton structure of the initial-state hadrons, such as the parton distribution functions (PDFs) and the transverse momentum-dependent distributions (TMDs) of proton. In recent years, heavy quarkonium production is also studied in heavy ion collisions to probe the quark-gluon plasma (QGP). The heavy quark pair is first produced in hard scattering at the early stage of the collisions and then interacts with the QGP and hadronizes to the heavy quarkonium on its way out of the QGP. Therefore, a good understanding of the heavy quarkonium production mechanism could facilitate our understanding of all these QCD objects.

A lot of data for the quarkonium production in different high-energy collisions have been collected. Take as an example, the cross section of ( collision) has been measured by the Belle and BaBar collaborations, the cross section of ( collision) has been measured by DELPHI collaboration at LEP, the yield and polarization of production in (photoproduction) have been measured by and Zeus at HERA, and the yield and polarization of in hadroproduction ( or collision) have been measured by CDF at Tevatron, by PHENIX and STAR at RHIC, and by CMS, ATLAS, ALICE, and LHCb experiments at the LHC. In the meantime, lots of theoretical efforts have been made to explain these experimental measurements. In this article, we will briefly review the current status of the phenomenological study of quarkonium production. We first overview theoretical frameworks for describing the inclusive quarkonium production. Then, we review the current status of the comparison between theoretical predictions and experimental measurements. We put stress on charmonium production at hadron colliders, with brief overview of the quarkonium production in , , and collisions. We refer the readers to several relevant reviews [15] for detailed discussions of other topics in quarkonium physics.

2. Quarkonium Production Mechanism

Heavy quarkonium production is usually separated into two steps: (1) the production of a pair with definite spin and color state in a hard collision, which could be calculated perturbatively, and (2) hadronization of the pair into a physical heavy quarkonium at a momentum scale much less than the heavy quark mass , which is in principle nonperturbative. Different treatments of the nonperturbative transition from pair to the physical quarkonium lead to different theoretical models. In the following, we briefly describe some of the most widely used ones: the color evaporation model (CEM) [68], the color-singlet model (CSM) [911], the nonrelativistic QCD (NRQCD) factorization theory [12], the fragmentation function approach [1317], and the most recently proposed soft gluon factorization (SGF) approach [18, 19].

2.1. The Color Evaporation Model (CEM)

In the CEM [68], it is assumed that every produced pair evolves into a specific heavy quarkonium if its invariant mass is below the open-charm/bottom threshold. It is further assumed that the probability for the pair to evolve into a specific quarkonium state is given by a constant which is independent of momentum and process. Mathematically, the production cross section of is expressed in the CEM as

where is the open-charm/bottom threshold. For each heavy quarkonium state, the CEM in equation (1) has one free parameter . The CEM is intuitive, simple, and successful to explain production data. However, it has a very strong prediction that the production rate of any two different charmonium states depends on neither the process nor kinematic variables, which contradicts data from many experiments. For example, the ratio of production cross section of to that of in collisions clearly depends on their transverse momentum [20, 21]. To overcome these obstacles, an improved version of the model, the ICEM, was proposed [22], in which the momentum of the pair is assumed to be larger than the momentum of quarkonium by a factor of . One consequence is that the lower limit of the above integral is replaced by . It was shown that the ICEM can describe the charmonium yields as well as the ratio of over [22]. The ICEM was also combined with -factorization to describe quarkonium polarization [2328].

2.2. The Color-Singlet Model (CSM)

In the CSM, the pair that evolves into the quarkonium is assumed to have the same color, spin, and orbital-angular-momentum quantum numbers as the heavy quarkonium. Particularly, it must be in a color-singlet state. Under this assumption, the production cross section for each quarkonium state is related to the wave function (or its derivatives) of around the origin, which can be extracted from the decay process of or calculated from the potential model or lattice QCD. Therefore, the CSM effectively has no free parameters. At relatively low energies, the LO CSM predictions for quarkonium production agree with the experimental data, while at high energies, the LO CSM predictions have been shown to underestimate the experimental data of direct and production at TeV collisions [29] by more than an order of magnitude, which is known as the surplus puzzle. In the past decade, it was found that the NLO and NNLO corrections to the CSM are significantly larger than the LO contributions [3032]. Including these corrections relieves the inconsistency between the LO CSM prediction and the data. However, a full description of data is still difficult. Besides, given the very large corrections at NLO and NNLO, it is not clear that the perturbative expansion in is convergent. Moreover, in the case of -wave production and decay, the CSM is known to be incomplete because it suffers from uncanceled infrared divergences. The last point can be rigorously cured in a more general framework of NRQCD factorization theory which we will discuss below.

2.3. The NRQCD Factorization Approach

NRQCD is an effective theory of QCD and reproduces full QCD dynamics at momentum scales of order and smaller. In NRQCD, the production cross section of a heavy quarkonium is given by the factorization formula [12]:

Here, is the NRQCD factorization scale, which is the ultraviolet (UV) cutoff of the NRQCD effective theory, is the short-distance coefficient (SDC) which describes the production of a pair with quantum number in the hard scattering, and is the NRQCD long-distance matrix element (LDME) that describes the hadronization of the pair in state into the heavy quarkonium . The LDME is defined as the vacuum expectation value of a four-fermion operator in NRQCD, and each LDME has a known scaling behavior in powers of . Then, the sum over can be organized in powers of . Therefore, equation (2) is a double expansion of and . In practice, for a certain accuracy, one truncates the summation and keeps only a few LDMEs for each production. The predictive power of the NRQCD factorization approach relies on the convergence of this velocity expansion, as well as the universality of LDMEs.

As shown in equation (2), the NRQCD factorization contains contributions from both the color-singlet (CS) and the color-octet (CO) channels. If one sets the CO contributions to zero, one could recover the CSM for -wave heavy quarkonium production. Thanks to the CO contributions, NRQCD solves the infrared divergence problem encountered in the CSM [33]. Although there is no all-order proof of NRQCD factorization for quarkonium production yet, it is found that the factorization holds at least to next-to-next-to-leading order (NNLO) in if the LDMEs are modified to be gauge complete [3438].

2.4. The Fragmentation Function Approach

The SDCs in the NRQCD factorization formula in equation (2) suffer from large high-order corrections for heavy quarkonium produced at large transverse momentum , which is an important kinematic region in high-energy colliders. In this region, the high-order corrections of the SDCs receive huge power enhancement in terms of , as well as large logarithmic corrections in terms of . To overcome these problems, a new QCD factorization approach was proposed to describe the large heavy quarkonium production [1317]. In this approach, the cross section is first expanded by powers of . Both the leading-power (LP) term and next-to-leading-power (NLP) term of the expansion could be factorized systematically into parton-production cross sections convoluted with several universal FFs, i.e.,

where the first term on the right side gives the contribution of LP in and the second term gives the NLP contribution. The symbol represents the convolution of the light-cone momentum fraction . In the first term, is the semi-inclusive cross section for initial hadrons and to produce an on-shell parton . The FFs represents the possibility of finding in the hadronization products of parton . The NLP contribution is similar, with an intermediate heavy quark pair in state instead of a single parton .

Since the NLP term is suppressed by , it seems to be not important at large . However, it is natural to expect that a heavy quark pair is more likely to evolve into a heavy quarkonium , comparing to a single quark or gluon. Consequently, the double-parton FFs are more important than the single-parton ones. This nonpertubative enhancement could balance the perturbative suppression in the intermediate range. There are more NLP contributions from other intermediate double partons besides a heavy quark pair. They are not included in equation (3) because they do not have this nonperturbative enhancement.

The factorization formula in equation (3) is a double expansion of and . The factorization scale is chosen at the same order of so no large logarithm exists. The FFs with different are related by a closed set of evolution equations. To make a theoretical prediction, one still needs a set of input FFs at a certain scale . These input FFs are nonperturbative and, in principle, should be extracted from fitting experimental data. However, since , it is natural to use NRQCD factorization to further factorize these input FFs. By doing this, all unknown input FFs could be expressed in terms of NRQCD LDMEs with the perturbatively calculable coefficients

where and describe the perturbative evolution of a parton and a pair with quantum number into a pair in the state , respectively. Mathematically, the FF+NRQCD factorization formula in equations (3), (4) and (5) is a reorganization of terms in equation (2). Physically, the FF+NRQCD factorization method correctly includes the evolution of a heavy quark pair when the relative velocity in the quarkonium rest frame is not much smaller than 1 and the NRQCD.

During the past two decades, the FFs in equation (4) and (5) have been widely studied. The coefficients for all double-parton FFs to both -wave and -wave states are calculated up to in Refs. [3941]. The coefficients for all single-parton FFs are available up to [4252] (see [3941] for a summary and comparison). At order, the coefficients of gluon FFs to , , , and are calculated in Refs. [38, 44, 5360]. Recently, the heavy quark FFs to state are obtained in Refs. [61, 62].

With the input FFs and the evolution equations, FF+NRQCD factorization formalism provides a systematic reorganization of the cross section in terms of powers of and a systematic method for resumming the potentially large -type logarithms. It is expected to have a better convergence in the expansion than NRQCD.

2.5. The Soft Gluon Factorization Approach

As we will discuss later, the NRQCD factorization still encounters some difficulties in describing inclusive quarkonium production data. It is known long time ago that NRQCD has bad convergence in velocity expansion [63], which may be responsible for the phenomenological difficulties. The aim of SGF is to provide a framework with better convergence [18].

In the SGF approach, the differential cross section of the quarkonium production is factorized as

where is the momentum of the intermediate pair, is the momentum of , and is soft gluon distribution function (SGD), which describes the hadronization of the pair into physical quarkonium by emitting soft hadrons. To account for the effect of soft hadron emission, which are mainly soft gluons perturbatively, the momentum of the intermediate state is kept different from the observed quarkonium momentum , which is different from the treatment in NRQCD. The SGD is defined by QCD fields in a small loop momentum region. With an explicit definition of the small region and taking advantage of equations of motion, the SGF is shown to be equivalent to the NRQCD factorization [19]. Nevertheless, compared with NRQCD, the SGF resums a series of relativistic corrections originating from kinematic effects, which results in a better convergence in velocity expansion. It is expected that the SGF approach may provide a better description of experimental data.

The first phenomenological application of the SGF approach was carried out in Ref. [64] for exclusive quarkonium production processes. It was shown there that, for production at factories, the SGF provides the best description of experimental data among all existing theoretical calculations. Recently, a quarkonium fragmentation function in SGF has been calculated to NLO in Ref. [65], which is the first step to apply SGF to inclusive quarkonium production processes. With explicit NLO calculation, it was demonstrated that the SGF is valid at the NLO level. Phenomenological application for inclusive quarkonium production in the SGF approach is still missing.

3. Quarkonium Production in Collisions

3.1. High Heavy Quarkonium Production

Based on the NRQCD factorization framework, the heavy quarkonium production in collisions has been widely studied. In the large region, the differential cross section of the quarkonium production can be factorized as

where ’s are the parton distribution functions (PDFs) for the partons in the initial colliding protons, is the collinear factorization scale, and is the partonic differential cross section.

At LO in , the partonic differential cross sections of are scaled as , , and , respectively, which are more important than the CS contribution , scaled as . By including the CO contributions, the surplus puzzle in CSM is solved naturally [66]. But NRQCD factorization encounters difficulties with charmonium polarizations. Since the dominant contribution is from the transverse channel, LO NRQCD predicts that and produced at hadron colliders are mainly transversely polarized [42, 67, 68]. On the contrary, experimental measurements at Tevatron and LHC find that these states are almost produced unpolarized [6971]. In addition, the LO calculation in NRQCD is difficult to explain the observed cross section ratio of the -wave charmonia at Tevatron [72], as the channel dominance predicts the ratio to be by spin counting [73, 74], which is much larger than the measured value of .

In the past decade, the perturbative SDCs for quarkonium production cross sections and polarizations have been calculated to NLO by three groups [7581]. At this order, the and channels are scaling, so their contributions are also important at large . Even though the NLO SDCs obtained by the three groups agree, they give very different predictions for polarization, due to the different methods used in fitting the CO LDMEs.

In Refs. [7577], it was found that at high , the SDC of the channel can be decomposed into a linear combination of the other two CO channels:

where for the Tevatron and for the LHC. Therefore, two linearly combined NRQCD LDMEs are introduced to fit the data:

where is or . Using the Tevatron data [82, 83] with GeV, the LDMEs are extracted as

With these fitted LDMEs, the theoretical predictions for prompt production are generally consistent with experimental data from the Tevatron and the LHC [8284] (see Figure 1). As shown in Figure 1(b), the polarization is roughly consistent with the CDF data, in which the is approximately unpolarized. Detailed analysis shows that the transversely polarized contributions from and channels cancel each other. Thus, a possible mechanism is that the unpolarized channel dominates, which leads to an unpolarized production [77, 85, 86]. For the production, such cancellation is weak [87], and its polarization is still hard to explain.

In Ref. [79], the authors determine the CO LDMEs by a global fit of a number of measurements. Using prompt production data in () [8894], ( GeV) [9597], [98], and [99] collisions, they determine all three CO LDMEs. In particular, they obtained a negative value for , and thus, the contributions from and channels add and enhance the transverse polarization at large . With these inputs, it was found that the is transversely polarized at large as shown in Figure 2.

In Ref. [80], the LDME determination is based on the yield data from CDF [89] and LHCb [94]. Similar to the method in Refs. [75, 77, 87], the data with are not considered in the fitting. It is found that both and are negative. As the two LDMEs have the same sign, the cancellation between transversely polarized contributions still occurs, which explains the unpolarized produced at large .

Due to the approximate heavy quark spin symmetry (HQSS) of NRQCD [12], the production of is closely related to the production of by the following relations:

Then, the measurement of can provide a further test of the LDMEs. The first prompt hadroproduction cross section was measured in 2014 by the LHCb collaboration with the centre-of-mass energies at  TeV and  TeV [100]. With the three sets of the LDMEs above, the authors of Ref. [101] found that theoretical calculations overshoot experimental measurement, as shown in Figure 3. Detailed analysis shows that the LHCb data are almost saturated by the contribution from the CS channel. This gives a constraint on . By assuming that the data are completely contributed from the channel, the authors of Ref. [102] obtain an upper bound for :

This bound is consistent with a previous study by the same authors on the yield and polarization [87], so they argue that the prompt production of and can be understood in the same theoretical framework. The authors of Ref. [103] fit both the CO and CS LDMEs, and their predictions are also compatible with data. Recently, the LHCb collaboration reports their measurement of the hadroproduction cross section at  TeV [104], which is consistent with the previous data.

Within the same framework as that for the production, the study of production at NLO was carried out by two groups [81, 105, 106]. Their results show a slightly transverse polarization at large , consistent with the measurements by CMS [107] within experimental uncertainties.

For -wave quarkonia, the first complete NLO study of production was performed in Ref. [108] in 2010. At NLO, the CS channels are scaled as , which give large contributions at high . They also find that decreases slower than , so the measured ratio of at the Tevatron [72] can be naturally explained (see Figure 4). In 2016, the authors of Ref. [109] performed a global analysis of the existing data on hadroproduction from the Tevatron [72] and the LHC [110113]. In the meantime, the polarization of the was also predicted [114116], which indicated that the polarization parameter of the should be positive while that of the should be negative. The experimental measurement was not available until 2019 by the CMS collaboration [117]. Current experimental data seem to be consistent with the NLO results.

As mentioned in the previous section, in the high region, the partonic differential cross sections in the NRQCD factorization formula in equation (7) contain large logarithms like , which could ruin the convergence of perturbative expansion. Thus, resummations of these large logarithmic terms are necessary. This can be done by using FF+NRQCD factorization approach. Applying the FF+NRQCD factorization approach with LP approximation was carried out in Refs. [85, 118] for the hadroproduction of , , and at high , and good agreement with the measurements is obtained. It was also found that contributions from and channels should almost cancel with each other, so the produced is almost unpolarized, which confirms the conclusion in Refs. [77, 87]. This implies that the qualitative results in the NLO NRQCD calculations are not changed by LP resummation. It is an interesting question whether the NLP resummation could change this conclusion.

3.2. Low Heavy Quarkonium Production

At the low ( is the quarkonium mass) region, the collinear factorization formalism in equation (7) is no longer applicable. At this regime, large (small , , where is the collider center of mass energy) contribution arises at higher orders that may not be fully accounted for in the collinear factorization framework. Another source of contribution is from the higher-twist multiparton matrix elements that are large at low . Both of these contributions can be computed systematically in the Color Glass Condensate (CGC) effective field theory [124, 125]. By combining the CGC and NRQCD formalisms, a new factorization framework for quarkonium production was proposed [123, 126], in which the SDCs in equation (2) are given by for the color-octet channels and

for the color-singlet channels. In these expressions, () is the transverse momentum (rapidity) of the produced heavy quarkonium, () is the rapidity of gluons coming from dilute proton (dense proton), denotes the fundamental dipole amplitude, and are the hard parts, is an unintegrated gluon distribution inside the proton, and is the effective transverse area of the proton. Such a CGC+NRQCD framework provides a good description of the yield [123] and polarization at low [127] in collisions, shown in Figures 5 and 6. Interestingly, the CGC+NRQCD result at small merges smoothly with the NLO NRQCD result at intermediate and large , providing an unified description for quarkonium production in the full region. Note that this framework is even more useful for quarkonium production in proton-nucleus collisions [128131].

4. Quarkonium Production in Annihilation at Factories

The quarkonium production in annihilation at factories is another important process to test the NRQCD factorization. For exclusive double charmonium production such as , the first theoretical calculation at LO in both and [135137] gives a production cross section about  fb. It is much smaller than the measured value, which is  fb by Belle [138] and  fb by BaBar [139], where is the branching fraction for the decaying into at least two charged tracks. Later, the authors of Refs. [140, 141] show that the NLO QCD correction can substantially enhance the cross section with a factor (the ratio of NLO to LO) of about . Meanwhile, the relative correction is also found to be significant [135, 142, 143]. Including both the and corrections may resolve the large discrepancy between theory and experiment. Recently, the NNLO QCD correction to has been completed [144], which gives a state-of-the-art calculation consistent with the BaBar measurement.

For inclusive production, the first measurements of the cross section are released by the BaBar [145] and Belle [138, 146] collaborations. It was found by Belle [138] that the cross section  fb is about a factor of larger than the LO NRQCD factorization predictions including both the CS [147151] and CO [151] contributions. Later, the Belle collaboration reported an updated measurement: [99], which is smaller than the previous one. But it is still much larger than the LO NRQCD calculation. The large gap between experiment and theoretical calculations is reduced by including the NLO QCD correction, which substantially enhances the cross section with a factor of about [152].

To explain measured by Belle [99], the NLO QCD correction to the CS channel is calculated in Refs. [153, 154], which increases the LO result by about . As shown in Figure 7, the NLO CS contribution saturates the Belle measurement. The relativistic correction [155, 156] and the QED initial state radiation (ISR) effect [157] have also been considered, which enhance the LO cross section by a factor of and , respectively. Surprisingly, including all these corrections leads to the CS contribution somewhat above the measurement by Belle, leaving little or no room for the contribution of the CO channel . This provides an upper limit of the CO LDMEs. By assuming a vanishing CS contribution and including the NLO QCD correction to , the authors of Ref. [158] obtain

which is much smaller than the value of the CO matrix element extracted from hadron colliders, i.e., [75, 76, 87], bringing the universality of LDMEs into question.

The universality problem and the polarization puzzle discussed in the previous section are two outstanding problems in NRQCD factorization. A possible solution is to resum high-order relativistic corrections, since for charmonium, is not a small number. This is exactly the motivation of the SGF approach. A detailed and comprehensive phenomenological study for inclusive quarkonium production in the SGF framework could help to understand these problems. In the meantime, more experimental data with smaller uncertainties are also indispensable.

5. Quarkonium Production in and Collisions

In the photoproduction of charmonia in collisions at HERA, a quasi-real photon emitted from the incoming electron interacts with a parton from the proton and produces a pair that evolves into a charmonium state. There are two types of processes contributing to the photoproduction cross sections. The first is the direct photoproduction, in which the virtual photon interacts with the parton electromagnetically. The second is the resolved photoproduction, in which the virtual photon emits a parton, which then interacts with the parton . Combining collinear factorization and NRQCD factorization, the inclusive photoproduction cross section can be written in the form of [159, 160]

where is the photon flux function, is either or the PDF of parton in the resolved photon, is the PDF of the proton, and is the partonic cross section.

The NLO QCD correction to the partonic cross sections for the direct production was first performed by Krämer in 1995 [161]. In 2009, complete NLO calculation was obtained in Refs. [159, 162, 163]. As shown in Figure 8, the yield data from HERA are roughly consistent with the global fit at NLO performed in Refs. [78, 160].

In addition to the unpolarized yield, the polarization of photoproduction has also been studied [162, 163, 165]. Figure 9 shows the LO and NLO NRQCD results of polarization observables as a function of by Ref. [165] comparing two HERA datasets. It is found that NRQCD predicts the produced at large to be approximately unpolarized, both at LO and NLO, which is confirmed by the data (Figure 9(a)). However, the ZEUS measurement (Figure 9(b)) exhibits a tendency towards transverse polarization. A possible explanation is that diffractively produced vector mesons prefer to be strongly transversely polarized in the endpoint region [165].

The inclusive cross section for production in colliders via fusion has also been measured by the DELPHI collaboration [98]. Similar to the photoproduction, both direct photon and resolved photon contribute to the cross section. The first complete NLO NRQCD computation including the resolved contributions was studied in Ref. [160]. However, by using the LDMEs obtained from the global fit, the predicted cross sections are several times below the DELPHI data, as shown in Figure 10. This is caused by a cancellation between the and contributions owing to the negative value of obtained by the global fit [2]. Thus, a global analysis of the world data is still challenging.

6. Summary

In this article, we have reviewed some theoretical methods to describe heavy quarkonium production, including the color evaporation model, color-singlet model, NRQCD factorization, fragmentation function approach, and soft gluon factorization. We then emphasize the current status of the phenomenological study of charmonium production, mainly in the NRQCD factorization framework. We concentrate on the comparison between theoretical predictions and experimental measurements for the charmonium production in four important processes: collision, annihilation, collision, and collision. After the NLO contributions are taken into account, the NRQCD factorization can give a qualitatively correct description for quarkonium production. In particular, by combining NRQCD factorization with CGC effective theory, a unified description for quarkonium hadroproduction in the full region is obtained. However, LDMEs determined from different choices of datasets can disagree with one another, and none of them are able to give a global description of all important observables, such as the total yield, momentum differential yield, and polarization. The polarization puzzle and the universality problem are two outstanding problems in NRQCD factorization, which may be caused by the bad convergence of velocity expansion. Hopefully, these difficulties could be resolved or relieved in the SGF framework with well-controlled relativistic corrections.

Disclosure

An ArXiv has previously been published [166].

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The work is supported by the National Natural Science Foundation of China (Grant Nos. 11875071 and 11975029), the National Key Research and Development Program of China under Contract No. 2020YFA0406400, and the Qilu Youth Scholar Funding of Shandong University.