Inclusive Hadroproduction of $P$-wave Heavy Quarkonia in pNRQCD

We compute the color singlet and color octet NRQCD long-distance matrix elements for inclusive production of $P$-wave quarkonia in the framework of pNRQCD. In this way, the color octet NRQCD long-distance matrix element can be determined without relying on measured cross section data, which has not been possible so far. We obtain inclusive cross sections of $\chi_{cJ}$ and $\chi_{bJ}$ at the LHC, which are in good agreement with data. In principle, the formalism developed in this work can be applied to all inclusive production processes of heavy quarkonia.

The mechanism underlying heavy quarkonium production is a key to understanding the dynamics of strongly coupled systems [1][2][3][4][5]. Quarkonium production is extensively studied in experiments at particle colliders like LHC, SuperKEKB, BEPC II, and RHIC, and will continue to be an important subject in future colliders such as the planned Electron-Ion Collider. Quarkonium production has a large impact on studies of the QCD phase diagrams and early universe, as the production in protonproton collisions is the bottom line to which quarkonium suppression in heavy ion collisions is compared [6]. Moreover, from the theoretical point of view, quarkonium production processes have exquisite theoretical issues pinning down on factorization in strongly coupled theories, definition and calculation of nonperturbative matrix elements, and resummation of logarithms of large ratios of scales [7][8][9][10][11].
The typical hierarchy of energy scales that characterizes heavy quarkonium is m mv mv 2 , where m is the heavy quark mass and v 1 the relative velocity of the quark in the bound state. This hierarchy of energy scales may be exploited to construct a hierarchy of effective field theories. Nonrelativistic QCD (NRQCD) [7,12] follows from QCD by integrating out modes associated with the energy scale m from Green's functions describing a heavy quark and a heavy antiquark near threshold. The matching to NRQCD can be done perturbatively, since m is larger than the typical hadronic scale Λ QCD . Potential NRQCD (pNRQCD) [13][14][15] follows from NRQCD by integrating out gluons of energy or momentum of order mv. The matching to pNRQCD may need to rely on nonperturbative methods if the momentum scale, mv, is comparable to Λ QCD . While NRQCD had great success in heavy quarkonium phenomenology, a satisfactory description of inclusive production processes from first principles is still beyond reach. Much of the difficulty stems from our limited knowledge of the NRQCD long-distance matrix elements (LDMEs), which describe the nonperturbative evolution of the heavy quark Q and antiquarkQ into a quarkonium. First-principles determinations have not been possible, even approximately, for a class of important LDMEs that are associated with the QQ in a color octet state. On the other hand, phenomenological determinations of the unknown LDMEs based on different choices of observables have led to inconsistent sets of LDMEs, which have resulted in contradicting predictions, in particular, leaving open the long-standing problem of the polarization of quarkonium produced in hadron colliders [16]. It would be of enormous impact to be able to compute the unknown LDMEs from first principles.
Potential NRQCD has been successfully applied to annihilation and exclusive electromagnetic production processes of heavy quarkonia [17][18][19]. It has been anticipated that pNRQCD could also be used to describe inclusive production processes. In this Letter, we apply for the first time pNRQCD to this kind of processes by computing the NRQCD LDMEs that appear in the inclusive production cross section of P -wave quarkonia. Specifically, we consider production cross sections of χ QJ (Q = c or b, J = 0, 1, and 2) at leading order in v.
The cross section is given in the NRQCD factorization formalism at leading order in v by [7] Here, we use spectroscopic notation for the angular momentum state of the QQ, while the superscripts 1 and 8 denote the color state of the QQ: color singlet (CS) and color octet (CO), respectively. The quantities σ QQ( 3 P [1] and σ QQ( 3 S [8] 1 ) are the perturbatively calculable shortdistance coefficients (SDCs). We have used the heavyquark spin symmetry to reduce the χ QJ LDMEs into LDMEs of χ Q0 , which are defined by where |Ω is the QCD vacuum, T a are SU(3) generators, σ i Pauli matrices, and ψ and χ Pauli spinors that annihilate and create a heavy quark and antiquark, respectively. The covariant derivative ← → D is defined by and A the gluon field. The operator P Q(P ) projects onto a state consisting of a quarkonium Q with momentum P . The path-ordered Wilson line along the spacetime direction , defined by where A adj is the gluon field in the adjoint representation, ensures the gauge invariance of the CO LDME [8]. The direction is arbitrary.
We aim at expressing the LDMEs (2) in pNRQCD. We work in the strong coupling regime, where Λ QCD mv 2 . This condition is fulfilled by non-Coulombic, strongly coupled quarkonia, such as possibly the χ QJ . In order to compute the LDMEs in strongly coupled pNRQCD, we expand the NRQCD Hamiltonian in inverse powers of the heavy quark mass, and compute at each order in 1/m in quantummechanical perturbation theory (QMPT) its eigenstates [20,21]: Here x 1 and x 2 are the positions of the quark and antiquark, respectively. Only for the eigenstates in the static limit, |n; x 1 , x 2 (0) , the quark and antiquark positions are conserved. The corresponding eigenvalues, E n (x 1 , x 2 ), are the energies of a static quark-antiquark pair located in x 1 and x 2 in the presence of light degrees of freedom in the ground state (n = 0) or in some excited state (n > 0). The set of quantum numbers denoted with n labels the excitations of the light degrees of freedom. We also define the states |n; x 1 , x 2 , which encode the light degrees of freedom content of the states |n; x 1 , x 2 , through the relation |n; The computation of the LDMEs in Eqs. (2) requires the knowledge of the matrix elements of P Q(P =0) on the states |n; x 1 , x 2 , which follows from general properties of the operator. We first note that P Q(P ) commutes with the NRQCD Hamiltonian. This is because the decay modes that involve momentum transfers at the scale m have been integrated out, and hadronic and electromagnetic transitions between quarkonium states can be ignored as they occur at a much larger time scale than the hadronization of Q andQ into quarkonium. Hence, the quarkonium state Q is conserved, and so is P Q(P ) . We write an eigenstate of H NRQCD and P Q(P ) as where φ Q(n,P ) (x 1 , x 2 ) is a suitable function of x 1 , x 2 , and P . We take |Q(n, P ) to be nonrelativistically normalized. For n = 0, the state is just a quarkonium state |Q(P ) [18]. For n = 0, the state includes excitations of the light degrees of freedom.
For the state |n; x 1 , x 2 to have a nonvanishing overlap with P Q(P ) , its static component |n; Tr being the trace over color and N c the number of colors [22]. This condition follows from requiring that, in the absence of gluonic excitations carried by 1/m corrections, the quark-antiquark fields at the same point are in a color-singlet configuration, which, in turn, is a necessary requirement for the state |n; x 1 , x 2 to overlap with a state containing a quarkonium. We denote with S the subset of eigenstates of the NRQCD Hamiltonian, |n; x 1 , x 2 , that fulfill the condition. Finally, the operator P Q(P ) can be written as In strongly coupled pNRQCD [14,20,21], color-singlet heavy quark-antiquark pairs in the presence of light degrees of freedom in a state n are described by the field S n and the Hamiltonian: is determined by matching order by order in 1/m to n; x 1 , x 2 |H NRQCD |n; x 1 , x 2 . In particular, at leading order in v, h n is the sum of a kinetic energy operator, p 2 /m, where p = −i(∇ 1 − ∇ 2 )/2 is the relative momentum of the heavy quark-antiquark pair, and a static potential that matches the eigenvalue E NRQCD . Since the eigenvalues of h n have to be equal to Q(n, P )|H NRQCD |Q(n, P ) , as pNRQCD and NRQCD describe the same physical spectrum, we can identify the functions φ Q(n,P ) (x 1 , x 2 ) with the eigenfunctions of h n . At leading order in v, it holds (8) where exp[−iP · (x 1 + x 2 )/2] is a plane wave describing the center of mass motion and φ We are now in the position to express the production LDMEs in terms of gluon field correlators and the wavefunctions at the origin. The procedure is similar to the one developed to compute the annihilation rates of heavy quarkonia in pNRQCD in Refs. [17][18][19] and consists of the following steps: (i) replace in the LDMEs the projector P Q(P =0) with the expressions (6) and (5); (ii) using QMPT, and in particular Eqs. (3) and (4), express the LDMEs in terms of |n; x 1 , x 2 (0) and E n (x 1 , x 2 ); (iii) make explicit the heavy quark and antiquark field content of the states |n; x 1 , x 2 (0) and eliminate the fields by using Wick's theorem; one makes use at this point of the fact that the states in P Q(P =0) belong to the set S, which constrains their color structure; (iv) rewrite the sum of the matrix elements of the gluon fields on the states |Ω and |n; x 1 , x 2 (0) (evaluated at x 1 − x 2 = 0) in terms of gluon field correlators; (v) identify φ Q(n,P ) (x 1 , x 2 ) or derivatives of them (evaluated at x 1 − x 2 = 0) with the wavefunctions of the pNRQCD Hamiltonian. The leading order wavefunctions can be computed by solving the corresponding Schrödinger equations, once the static potentials have been determined from the static energies E  Q (x 1 , x 2 ) that correspond to the n = 0 case. For n = 0, the static potential in h 0 can be extracted from the vacuum expectation value (VEV) of a static Wilson loop [14,20]. Similarly, for n = 0 and n ∈ S, the static potential in h n can be extracted from the VEV of a static Wilson loop in the presence of some additional, disconnected gluon fields. To our knowledge, there are no lattice data available for the n = 0 static potentials. However, we expect that the disconnected gluon fields mostly provide a constant shift to the potentials, for instance in the form of a glueball mass, but do not significantly affect their slopes. This is also supported by large N c considerations. In the large N c limit, the VEV of a Wilson loop with additional disconnected gluon fields factorizes into the VEV of the Wilson loop times the VEV of the additional gluon fields up to corrections of order 1/N 2 c [23,24]. If the slopes of the static potentials are the same for all n in the large N c limit, then in that limit the wavefunctions φ (0) Q(n) (x 1 , x 2 ) are independent of n. Hence, we will approximate the wavefunctions φ  1 ) in strongly coupled pNRQCD. Furthermore, in the case of the CO LDME, we approximate φ where R (0) χ Q0 (r) is the radial wavefunction of χ Q0 at leading order in the velocity expansion (R (0) χ Q0 (r) stands for its derivative). This reproduces the result obtained in the vacuum-saturation approximation in Ref. [7].
The CO LDME O χ Q0 ( 3 S [8] 1 ) vanishes at leading order in QMPT. Nonvanishing contributions come from nextto-leading order in QMPT: E a,i (t) being a chromoelectric field component computed at the time t and at the space location 0, and Φ 0 (t, t ) = P exp[−ig t t dτ A adj 0 (τ, 0)] a Schwinger line. Note that E is a purely gluonic quantity that does not depend on the heavy quark flavor.
The expression for the CO LDME given in Eq. (10) is very similar to the pNRQCD expression for the CO LDME appearing at leading order in v in the decay of χ QJ into light hadrons [17,19]. The only difference is that there the correlator E is replaced by the correlator The two correlators would be the same if we could neglect the contributions from the strings. At one loop, they have the same logarithmic dependence on the renormalization scale Λ and satisfy the same evolution equation. From this [17], it follows that the one-loop evolution equation where C F = (N 2 c − 1)/(2N c ). Equation (13) agrees with the evolution equation derived from a perturbative calculation in NRQCD [7]. The agreement is a one-loop consistency check of Eq. (10). At two loops however the identification of E with E 3 may not hold [8][9][10][11].
Equation (10) is our result for the CO LDME. The result allows a first-principles determination of the CO LDME, once E is known. The correlator E may be computed in lattice QCD or it can be obtained from processes involving heavy quarkonia.
We now compute the inclusive production cross sections of χ cJ and χ bJ (nP ) from proton-proton collisions at the LHC based on our results for the LDMEs in Eqs. (9) and (10). We use the value |R (0) χc0 (0)| 2 = 0.057 GeV 5 , which we obtain by comparing the measured two-photon decay rates of χ c0 and χ c2 in Ref. [25] with the pNRQCD expressions at leading order in v and at next-to-leading order (NLO) in α s [19]. Because two-photon decay rates of χ bJ have not been measured yet, we take for |R (0) χ b0 (nP ) (0)| 2 the averages of the values listed in Ref. [19], which are obtained from the potential-model calculations in Refs. [26][27][28][29]. We take |R χ Q0 (0)|, E, and m. We include the uncertainties coming from E and from unaccounted relativistic corrections of relative order v 2 , which we take to be 30% (10%) of the central values of the χ cJ (χ bJ ) cross sections. We neglect the uncertainty from the corrections of order 1/N 2 c , which is small compared to other uncertainties. We add the uncertainties in quadrature. 15 20 25 30 FIG. 1. Differential cross sections for χc1 and χc2 at the LHC ( √ s = 7 TeV, |y| < 0.75), compared with the ATLAS measurements [30].
We fit E on χ cJ production data. In particular, we compute production cross sections of χ c1 and χ c2 by taking the SDCs in Ref. [31] that are accurate at NLO in α s and include resummed leading logarithms of p T /m for √ s = 7 TeV in the rapidity range |y| < 0.75. Here, p T is the transverse momentum of the χ cJ . We present the differential cross sections B × dσ/dp T , where B = Br(χ cJ → J/ψγ) × Br(J/ψ → µ + µ − ) is taken from Ref. [32]. Our results for the differential cross sections are shown in Fig. 1 against ATLAS data [30]. The resulting fitted value of E in the MS scheme is The fit quality is good. We note that the value (14) is close to E 3 (Λ = 1.5 GeV) = 2.73 +0.94 −0.65 obtained from χ cJ decay data [19]. This, together with the 1 loop running agreement, is suggestive that for the considered P -wave CO LDME at leading order in v an approximate crossing relation may hold.
In order to compare with data also in the bottomonium case, we compute the feeddown fractions R Υ(n S) at the LHC ( √ s = 7 TeV, 2 < y < 4.5), compared with the LHCb data [33].
which are defined as the fractions of Υ(n S) (n = 1, 2, and 3) originating from decays of χ b1 (nP ) and χ b2 (nP ) (n =1, 2, and 3). We take the correlator E extracted from the charmonium production data and run it to the bottom mass at one loop accuracy. We compute the SDCs at NLO in α s for √ s = 7 TeV in the rapidity range 2 < y < 4.5 by using the FDCHQHP package [34]. We compute R χ b (nP ) Υ(n S) by multiplying the χ bJ (nP ) cross section by the branching ratio Br(χ bJ (nP ) → Υ(n S) + γ), summing over J = 1 and 2, and dividing by the inclusive Υ(n S) cross section. For the branching ratios of the 1P and 2P states, we take the measured values from Ref. [32], while, for the branching ratios of the 3P states, we take the theoretical predictions from Ref. [35] because there are no data available. We present our results as functions of the transverse momentum p Υ(n S) T of the Υ(n S), which we compute by multiplying the transverse momentum of the χ bJ (nP ) by m Υ(n S) /m χ bJ (nP ) . We take the bottomonium masses from Ref. [32]. For the denominator, we take the Υ(n S) LDMEs determined in Ref. [35], which give good descriptions of the measured cross sections at the LHC for p T > 15 GeV. Our results for R χ b (nP ) Υ(n S) are shown in Fig. 2 against the LHCb data at √ s = 7 TeV [33]. Our determinations are in good agreement with the data, with the possible exception of the χ b (1P ) case, which, however, may not be well described as a strongly coupled quarkonium (see e.g. [36]). Our determinations have been made possible by Eq. (10), with E fixed on χ cJ production data, while previous determinations, including the ones in Ref. [35], relied on using measured bottomonium cross section data to determine the χ bJ (nP ) LDMEs for each n.
The pNRQCD approach for inclusive production of Pwave heavy quarkonia that we have developed in this Letter provides expressions for the color-singlet and coloroctet NRQCD LDMEs in terms of quarkonium wavefunctions and universal gluonic correlators, which may be determined from data and/or lattice QCD calculations. This brings in a reduction in the number of nonperturbative unknowns and a substantial enhancement in the predictive power of the nonrelativistic effective field theory approach to inclusive heavy quarkonium production processes. Our results have made possible for the first time to determine the χ bJ (nP ) cross sections from first principles without fitting color-octet LDMEs on bottomonium production data. Our results for the χ cJ cross sections and χ b (nP ) feeddown fractions at the LHC have been summarized in Figs. 1 and 2. They are in good agreement with measurements. We note that the formalism developed in this work may be applied to all heavy quarkonium production processes and, in particular, to the production of the J/ψ or the Υ(nS) states [37]. It may be also possible to extend the formalism to inclusive production of heavy hybrids and mixed states of hybrids and quarkonia.
We are grateful to Geoffrey Bodwin for pointing out an inconsistency in an early version of this work.