Inclusive Production of Heavy Quarkonia in pNRQCD

We develop a formalism for computing inclusive production cross sections of heavy quarkonia based on the nonrelativistic QCD and the potential nonrelativistic QCD effective field theories. Our formalism applies to strongly coupled quarkonia, which include excited charmonium and bottomonium states. Analogously to heavy quarkonium decay processes, we express nonrelativistic QCD long-distance matrix elements in terms of quarkonium wavefunctions at the origin and universal gluonic correlators. Our expressions for the long-distance matrix elements are valid up to corrections of order $1/N_c^2$. These expressions enhance the predictive power of the nonrelativistic effective field theory approach to inclusive production processes by reducing the number of nonperturbative unknowns, and make possible first-principle determinations of long-distance matrix elements once the gluonic correlators are known. Based on this formalism, we compute the production cross sections of $P$-wave charmonia and bottomonia at the LHC, and find good agreement with measurements.


Introduction
Understanding the mechanism of inclusive heavy quarkonium production is one of the most challenging problems in QCD [1][2][3][4]. Heavy quarkonium production processes provide probes of the interplay between perturbative and nonperturbative aspects of QCD, and are considered key in understanding the hot and dense quark-gluon plasma. For decades, much theoretical effort has been made using the nonrelativistic QCD (NRQCD) factorization formalism [5]. In this formalism, inclusive production rates of a heavy quarkonium are factorized into products of short-distance coefficients, which encode the perturbative physics at the scale of the heavy quark mass and above, and long-distance matrix elements (LDMEs) that depend on the nonperturbative nature of the quarkonium state. While much progress has been made in computing the short-distance coefficients in perturbative QCD, it remains unknown how to compute from first principles a wide class of LDMEs. Therefore, in most phenomenological studies, the long-distance matrix elements have been obtained from fits to cross section data. This approach has resulted in different sets of LDME determinations depending on the choice of data that are used in the fit, which can disagree with one another [6]. These inconsistent sets of LDMEs lead to contradicting predictions, and none of the determinations is able to give a comprehensive description of the most important observables related to heavy quarkonium production at a satisfactory level [6,7]. It is therefore desirable and to some extent even necessary to have first-principle based constraints, or even computations, of the NRQCD LDMEs, in order to enhance significantly our understanding of the heavy quarkonium production mechanism.
Nonrelativistic effective field theories such as NRQCD exploit the hierarchy of energy scales that appear in processes involving heavy quarkonia, which are the heavy quark mass m, the typical relative momentum mv of the quark and antiquark, and the typical binding energy mv 2 [8]. Here, v 1 is the relative velocity of the heavy quark and the antiquark inside a heavy quarkonium. 1 NRQCD is obtained by integrating out modes associated with energy scales of order m and higher [5,9]. These are encoded in the short-distance coefficients, while all contributions from the scales mv and mv 2 are contained in the NRQCD LDMEs. Because the heavy quark mass is much larger than Λ QCD , the short-distance coefficients can be computed in perturbative QCD, and are given by series in α s . The LDMEs that appear in the NRQCD factorization formulas for inclusive production cross sections of heavy quarkonia describe the evolution of a heavy quark and antiquark pair into a heavy quarkonium state. This should be contrasted with the LDMEs for decay rates of heavy quarkonia that correspond to probabilities of finding a heavy quark and antiquark pair inside a heavy quarkonium. In case of LDMEs that involve heavy quark and antiquark pairs in color-singlet states, the so-called vacuum-saturation approximation can be used to relate the production and decay LDMEs. Hence, the color-singlet production LDMEs have been obtained from lattice QCD, potential-model calculations, or from measured decay rates. On the other hand, it is not known how to relate the production and decay LDMEs that involve heavy quark and antiquark pairs in color-octet states, and the color-octet production LDMEs have not been computed from first principles.
The effective field theory potential NRQCD (pNRQCD) [8,[10][11][12][13] is obtained by further integrating out modes associated with energy scales larger than mv 2 appearing in NRQCD. In the regime mv 2 Λ QCD , which is satisfied by non-Coulombic, strongly coupled quarkonia, NRQCD decay LDMEs can be factorized into products of quarkonium wavefunctions at the origin and universal gluonic correlators. The quarkonium wavefunctions are determined by solving the Schrödinger equation, which is the equation of motion of pNRQCD at leading order in v. The gluonic correlators are universal quantities that can be computed in lattice QCD. Strongly coupled potential NRQCD has been successfully applied to heavy quarkonium decay and exclusive electromagnetic production processes [14][15][16][17]. It has been anticipated that a similar formalism could be used to describe inclusive production processes of heavy quarkonia, by computing the production LDMEs in pNRQCD. Such a calculation would allow model-independent, first-principle based predictions of heavy quarkonium production cross sections.
Calculations of production LDMEs are also important in understanding the NRQCD factorization formalism for inclusive production of heavy quarkonia. The validity of the factorization formalism depends on the infrared finiteness of the short-distance coefficients to all orders in the expansion in powers of α s . An essential ingredient in proving such factorization is the determination of the infrared behavior of the production LDMEs. So far, the infrared properties of the production LDMEs have only been studied to two-loop accuracy [18][19][20][21][22]. It is possible that the pNRQCD expressions of the production LDMEs will simplify the investigation of their infrared properties and the verification of the NRQCD factorization formalism.
We have presented a first calculation of production LDMEs in strongly coupled pN-RQCD in ref. [23] for production of P -wave heavy quarkonia. In this paper, we describe in detail the formalism and clarify some technical details that may be useful in future studies. We work in the strong coupling regime mv 2 Λ QCD , and with heavy quarkonium states that are below the open flavor threshold. The calculation of the production LDMEs in this paper is valid up to corrections of relative order 1/N 2 c and v 2 . 2 Furthermore, we will present an extended set of phenomenological results.
The paper is organized as follows. In section 2, we set up the general formalism for the computation of the production LDMEs in pNRQCD. In section 3, we express in pNRQCD the production LDMEs for P -wave heavy quarkonium states. Phenomenological applications that include the computation of the cross sections for χ cJ (1P ) and χ bJ (nP ) states at the LHC can be found in section 4. We conclude in section 5.

LDMEs in pNRQCD
In the NRQCD factorization formalism, the inclusive production cross section of a quarkonium Q is given by [5] σ Q+X = N σ QQ(N ) Ω|O Q (N )|Ω . (2.1) Here, σ QQ(N ) are short-distance coefficients that correspond to the production cross section of a heavy quark-antiquark pair (QQ) in a spin and color state N , and |Ω is the QCD vacuum state. The NRQCD long-distance matrix element Ω|O Q (N )|Ω describes the evolution of the QQ in a state N into the quarkonium Q. The matrix elements have known scalings in v, which allows us to organize the sum over N in eq. (2.1) in powers of v. In general, the QQ can be in a color-singlet or in a color-octet state. For a color-singlet state, the operators O Q (N ) have the form where ψ and χ are Pauli spinor fields that annihilate and create a heavy quark and antiquark at spacetime position 0, respectively. The operator P Q(P ) projects onto a state consisting of a heavy quarkonium Q with momentum P ; it may be written as P Q(P ) = a † Q(P ) a Q(P ) , where a † Q(P ) is an operator that creates a quarkonium Q with momentum P . The quantities K N 2 In the special kinematical situation mv 2 ΛQCD mv, the ratio mv 2 /ΛQCD is larger than v and the induced corrections to the LDMEs may turn out to be parametrically larger than v 2 , which is their natural size when ΛQCD ∼ mv. In order not to complicate unnecessarily the error estimate, in the rest of the paper we will count the neglected corrections to the LDMEs in the v expansion as v 2 for all the kinematical situations that fulfill mv 2 ΛQCD. and K N are polynomials of Pauli matrices and covariant derivatives, and are proportional to the color identity matrix 1 c , so that the operators χ † K N ψ and ψ † K N χ are color singlets.
For a color-octet state, the operators O Q (N ) take the form [18][19][20] where T a is a color matrix in the fundamental representation. The operator Φ (x) is a Wilson line along the direction in the adjoint representation defined by where A adj is the gluon field in the adjoint representation and P stands for the path ordering of the color matrices. The Wilson lines Φ are necessary to ensure the gauge invariance of the vacuum expectation value of eq. (2.3), because in the absence of Φ , the gauge transformations of χ † K N T a ψ and ψ † K N T c χ do not commute with the operator P Q(P =0) . The direction is arbitrary. Hereafter we refer to the Φ in eq. (2.3) as gauge-completion Wilson lines. The NRQCD factorization conjecture is the statement that the short-distance coefficients in eq. (2.1) can be computed in perturbative QCD, and the NRQCD matrix elements are universal. The short-distance coefficients are computed in perturbative QCD by replacing the heavy quarkonium state Q in eq. (2.1) with a perturbative QQ state. For the short-distance coefficients to be perturbatively calculable, the infrared (IR) divergences that appear in the perturbative QCD calculation of eq. (2.1) must be reproduced by the NRQCD matrix elements to all orders in α s , so that the short-distance coefficients are IR finite. The outline for an all-orders proof of the IR finiteness of the short-distance coefficients has been given in refs. [19,24] for the case where the quarkonium Q is produced with a momentum that exceeds the mass of the heavy quarkonium m Q : for example, the proof applies to the case when the quarkonium is produced with a large transverse momentum p T m Q and the cross section is expanded to next-to-leading power in m 2 Q /p 2 T . On the other hand, the universality of the NRQCD matrix elements requires the color-octet matrix elements to be independent of the direction of the gauge-completion Wilson lines. A proof of the universality based on traditional methods of perturbative factorization has only been investigated to next-to-next-to-leading order in α s for the specific case where K N = K N = 1 [18][19][20][21][22].
In this paper, we aim at expressing the color-singlet and color-octet long distance matrix elements in pNRQCD. We work in the strong coupling regime, where mv 2 Λ QCD . This condition is fulfilled by non-Coulombic, strongly coupled quarkonia. The computation is based on expanding the NRQCD Hamiltonian in inverse powers of the heavy quark mass m, and on making use of quantum-mechanical perturbation theory to compute its eigenstates order by order in 1/m [12,13]. Because the energy scales that appear in NRQCD are mv, mv 2 , and Λ QCD , the expansion in powers of 1/m in the Hamiltonian leads to an expansion in powers of the dimensionless parameters v and Λ QCD /m in the observables. Explicitly, the static NRQCD Hamiltonian, H NRQCD , and its 1/m correction H (1) NRQCD are given by where E i a T a = E i = G i0 and B i a T a = B i = − ijk G jk /2 are the chromoelectric and chromomagnetic fields, respectively, G µν a T a = G µν is the gluon field strength tensor, q k are n f massless quark fields, and D = ∇ − igA is the gauge covariant derivative. The matrices σ i are the Pauli matrices, and c F is a short-distance coefficient, which is known to three-loop accuracy [25]. The physical states |phys are constrained by the Gauss law, which reads The space of states of NRQCD in the QQ sector is spanned by states |n; x 1 , x 2 containing a heavy quark located at x 1 and a heavy antiquark located at x 2 : where the states |n; x 1 , x 2 do not contain heavy particles. 3 We may take both |n; x 1 , x 2 and |n; x 1 , x 2 to satisfy orthonormality relations. Furthermore, the states |n; x 1 , x 2 are such that the NRQCD Hamiltonian is diagonal in n on them: Exploiting the expansion (2.5), the states |n; x 1 , x 2 can be computed order by order in 1/m: The states |n; NRQCD , whose eigenvalues are the static energies E (0) n (x 1 , x 2 ); the coordinates x 1 and x 2 are conserved in the static limit. The states |n; x 1 , x 2 (1) /m are the order-1/m corrections due to H NRQCD /m. The explicit expression of |n; x 1 , x 2 (1) in terms of |n; x 1 , x 2 (0) and E (0) n (x 1 , x 2 ) can be found in refs. [12,13]. 3 For further use we make explicit the color indices of the state: The number n labels the QCD static energies E (0) n (x 1 , x 2 ). Each static energy describes physically a distinct excitation of the static quark-antiquark pair due to the light degrees of freedom (gluons, light quarks). Lattice calculations of the QCD static spectrum suggest that each excitation is separated from the other by an energy gap of order mv or Λ QCD [26][27][28]. The n = 0 state corresponds to the ground state. For vanishing heavy quark-antiquark distance, r = x 1 − x 2 , the ground state reduces to the QCD vacuum: 12) where N c is the number of colors. This follows from the fact that the gluonic component of the Fock state describing a static heavy quark-antiquark pair in the ground state located at the origin is proportional to the color identity, the state being a color singlet, and is proportional to the QCD vacuum state, the state being the ground state. We compute the matrix elements Ω|O Q (N )|Ω by first evaluating the matrix elements of the operator P Q(P =0) on the states |n; x 1 , x 2 in pNRQCD, which we express in terms of quarkonium wavefunctions. Then, the remaining matrix elements of the heavy quark and gluon operators can be evaluated by using the methods first developed in refs. [14,15] for computing the decay matrix elements in pNRQCD.

Matching of P Q(P)
In this section, we express the matrix elements of the operator P Q(P ) = a † Q(P ) a Q(P ) between states |n; x 1 , x 2 in terms of quarkonium wavefunctions. Our computation is based on some general properties of the operator.
In NRQCD, processes that involve momentum transfers at the scale m are integrated out. Hence, production and decay of heavy quarkonium can occur within NRQCD only through hadronic and electromagnetic transitions. However, the time scales associated with these processes are much larger than the time scale associated with the hadronization of QQ into heavy quarkonium, which is typically of order 1/Λ QCD . Therefore, to a good approximation, we can neglect transition processes in describing the hadronization process within NRQCD. We conclude that the operator P Q(P ) , which counts the number of heavy quarkonia, is conserved. Since it commutes with the NRQCD Hamiltonian, P Q(P ) and H NRQCD can be diagonalized simultaneously. A simultaneous eigenstate of P Q(P ) and H NRQCD has the form (2.13) The states |n; x 1 , x 2 have been introduced above. They may transform as color singlets or color octets. Only color-singlet states project on Ω|χ † K N ψ and only color-octet states project on Ω|χ † K N T a ψΦ †ab (0). In |n; x 1 , x 2 we have dropped the color label because the spectrum, E n (x 1 , x 2 ; ∇ 1 , ∇ 2 ), does not depend on it. Color labels are implicit for coloroctet states. The functions φ Q(n,P ) (x 1 , x 2 ) are eigenfunctions of E n (x 1 , x 2 ; ∇ 1 , ∇ 2 ). They depend on the heavy quark and antiquark locations, on the quarkonium Q that moves with momentum P , and on the light degrees of freedom through the quantum number n. For n = 0 the state in eq. (2.13) is just a heavy quarkonium state [15]. We take the state |Q(n, P ) to be nonrelativistically normalized. Infinitely massive quark-antiquark pairs must be in a color-singlet state at the production point x 1 = x 2 in order to overlap with a state containing a quarkonium. This requirement limits the number of states that contribute to P Q(P ) . In order to contribute to P Q(P ) , states must have a static component |n; where Tr is the trace over the color indices. Property (2.14) will turn out to be crucial in evaluating the LDMEs, which, involving the local operators χ † K N ψ, ψ † K N χ, χ † K N T a ψ and ψ † K N T c χ, require the quark and antiquark to be located at the same point. We denote by S the subset of eigenstates that fulfill the requirement (2.14). 4 Then, the explicit expression of the projector P Q(P ) is It implies that n; x 1 , x 2 |P Q(P ) |k; x 1 , x 2 = 0 for n or k / ∈ S and for n = k, while n; x 1 , x 2 | P Q(P ) |n; x 1 , x 2 is completely determined by the function φ Q(n,P ) (x 1 , x 2 ) when n ∈ S. We will now show how to compute these functions at leading order in v.
Under the assumption mv 2 Λ QCD , the dynamics of the degrees of freedom scaling with the energy scale mv 2 is described by the Hamiltonian of strongly coupled pNRQCD. Strongly coupled pNRQCD follows from NRQCD by integrating out degrees of freedom carrying momentum or energy of order mv or Λ QCD [11][12][13]. It is assumed that its degrees of freedom consist of color-singlet states made of a heavy quark-antiquark pair and light degrees of freedom. If the effect of dynamical light quarks is neglected (the coupling of light mesons to quarkonia below threshold is a subleading effect, suppressed by powers of momenta of the light mesons, which is of order mv 2 ), then the dynamical degrees of freedom of pNRQCD consist of quarkonia and quarkonium exotica [14,15,17]. These degrees of freedom are represented by fields S n annihilating a color-singlet heavy quark-antiquark pair with light degrees of freedom in a state n. If we assume the simplest scenario that the energy levels for different n are separated by a gap of order Λ QCD , or that the mixing between fields with different n may be neglected, then we can isolate the dynamics of one single field S n . The strongly coupled pNRQCD Hamiltonian that describes this dynamics is The function h n (x 1 , is determined by matching to the NRQCD matrix element n; x 1 , x 2 |H NRQCD |n; x 1 , x 2 of eq. (2.10). At leading order in v, h n (x 1 , x 2 ; ∇ 1 , ∇ 2 ) is given by is the kinetic energy of the heavy quark-antiquark pair and V (0;n) (x 1 , x 2 ) is the static potential. The matching fixes the static potential to be the energy eigenvalue E (0) NRQCD . As a consequence of the matching, the functions φ Q(n,P ) (x 1 , x 2 ) are also eigenfunctions of h n .
For n = 0, V (0;n=0) is the QCD static potential, which can be determined from the vacuum expectation value of a static Wilson loop. In particular, it holds that where W r×T is a rectangular Wilson loop of spatial and temporal lengths r = |x 1 − x 2 | and T , respectively. Then, the function φ Q(n=0,P ) (x 1 , x 2 ) at leading order in v is given by where e iP ·(x 1 +x 2 )/2 is a plane wave encoding the center-of-mass motion, and φ For n ∈ S and n = 0, the static potential V (0;n) can be obtained from the vacuum expectation value of a static Wilson loop whose initial and final states select the excitation n of the static sources. While lattice QCD determinations of V (0;n) for n ∈ S and n = 0 are not available yet, we expect the disconnected gluon fields to produce mainly a constant shift to the potentials, possibly in the form of a glueball mass, and do not significantly affect the slopes. This is also supported by the large N c limit, where the vacuum expectation value of a Wilson loop with additional disconnected gluon fields factorizes into the vacuum expectation value of the Wilson loop times the vacuum expectation value of the additional gluon fields up to corrections of order 1/N 2 c [29,30]. If the slopes of the static potentials are the same for all n ∈ S, then the wavefunctions φ Q(n,P ) (x 1 , x 2 ) are given by eq. (2.19), and are independent of n. Hence, we take the approximation 20) making an error that is at most of order 1/N 2 c and of order v 2 . We can now compute the matrix elements of P Q(P =0) on the states |n; x 1 , x 2 . From eqs. (2.15), (2.13), and (2.20), we obtain This leads to the pNRQCD expression for the operator P Q(P =0) :

Matching of the LDMEs
Equation (2.22) implies that a LDME can be matched into the following pNRQCD expression The space integration is compensated by the denominator The contact terms in eqs. (2.24) and (2.25) can be computed by substituting the states |n; x 1 , x 2 with their quantum-mechanical expansion in powers of 1/m, see eq. (2.11), and by making explicit the heavy quark and antiquark fields using eq. (2.9). The heavy quark and antiquark fields are then removed by using Wick's theorem. The calculation of the contact terms is done in a similar way as the computation of the decay matrix elements in strongly coupled pNRQCD, except that the intermediate states are not just a quarkonium state, but include all states belonging to S and, in the calculation of the color-octet production matrix elements, the gauge-completion Wilson lines must be included. As a result, the contact term V O(N ) is proportional to the delta function δ (3) (r) or to derivatives of it. Furthermore, it can involve matrix elements of gluon fields that may be eventually expressed in terms of temporal correlators of gluon fields by using the techniques developed in refs. [14,15]. Finally, the LDMEs for inclusive heavy quarkonium production are computed from eq. (2.23). Because of the form of the contact terms, they turn out to depend on the quarkonium wavefunctions at the origin or its derivatives, and on vacuum expectation values of gluon fields.
The contact term for the color-octet matrix element in eq. (2.25) can in principle depend on the direction of the gauge-completion Wilson lines Φ . The dependence on the direction disappears in the LDMEs if the quarkonium Q has zero angular momentum, because the product of wavefunctions φ (0) in the integrand of eq. (2.23) is isotropic. The dependence also disappears in color-octet matrix elements for production of a quarkonium with nonzero angular momentum, if we sum over the quarkonium polarizations. The disappearance of the dependence is a necessary condition for the NRQCD factorization to hold [18,19]. Therefore, we can already state, on general grounds, that the calculation in pNRQCD of the LDMEs will support their universality in the case of unpolarized and polarization-summed production cross sections of heavy quarkonia.
In order to be consistent with the perturbative factorization of eq. (2.1), the righthand side of eq. (2.23), when computed in perturbative QCD, must have the same IR divergences as the NRQCD counterpart, which is obtained by computing in perturbative QCD the vacuum expectation values of the operators in eqs. (2.2) and (2.3). In section 3.2, we will confirm this agreement for the color-octet matrix elements that appear in the production cross sections of P -wave quarkonia.
3 Theory of inclusive production of P -wave quarkonia

P -wave LDMEs in pNRQCD
Based on the formalism developed in the previous section, we compute here the LDMEs that appear in production cross sections of P -wave heavy quarkonia, which include h Q and χ QJ , where Q = c or b and J = 0, 1, and 2. The NRQCD factorization formula for the inclusive production cross sections of P -wave quarkonia at leading order in v read The operators are where we have used the notation The parameter λ is the polarization of the quarkonium with nonzero angular momentum. We sum over all polarizations λ when computing polarization-summed cross sections, while we do not sum when computing polarized cross sections. In this section, we restrict to the case of polarization-summed cross sections.
Since we neglect transition processes between heavy quarkonium states in our treatment of the LDMEs, the inclusive production cross sections that we compute from the NRQCD factorization formula include only the "direct" production rates, where feeddown contributions that come from decays of higher quarkonium states are neglected. The feeddown contributions to inclusive quarkonium production cross sections can be included by adding direct production cross sections of higher quarkonium states, multiplied by the branching ratios into the measured quarkonium (see, for instance, section 4.4).
We first match the contact term for the operator 1 ). At leading order in the quantum-mechanical perturbation theory, it holds that where r = x 1 − x 2 and ∇ r = (∇ 1 − ∇ 2 )/2. In the last equality, we keep only terms that give nonzero contributions when inserted in eq. (2.23) to compute LDMEs for Pwave quarkonium states. Because the P -wave wavefunctions φ (0) 1 ) that contain two derivatives can make nonvanishing contributions to the right-hand side of eq. (2.23). We get terms proportional to derivatives only if the intermediate state is taken to be n = 0, an observation that follows from eq. (2.12), the canonical commutation relations and symmetry considerations [8,12]; the n = 0 state eventually reduces to the vacuum state so that Ω|D|Ω = ∇. From (3.3) we get which implies at leading order in v where the radial wavefunction , λ being the polarization of the 1 P 1 state and Y λ 1 (r) the spherical harmonics for P -wave states. The factor 3 in eq. (3.5) comes from the sum over the 3 polarizations of the h Q , while the factor 3/(2π) comes from the trace over the heavy quark spin times the . At leading order in v, the radial wavefunction is independent of the polarization λ. Equation (3.5) reproduces the well-known result obtained in the vacuum-saturation approximation [5].
In the following paragraphs, we match the contact term for the color-octet operator 0 ). At leading order in the quantum-mechanical perturbation theory, we have This expression leads to a vanishing contribution to the LDME once inserted in the righthand side of eq. (2.23). The reason is twofold: first, there are no derivatives acting on the P -wave wavefunctions and second, the heavy quark-antiquark pair in the state |n; behaves like a color singlet at r = 0 (see eq. (2.14)), which leads to a trace over an SU (3) generator after application of Wick's theorem.
In order to obtain a nonvanishing contribution to the octet LDME, we need to include 1/m corrections to the states |n; x 1 , x 2 appearing in the left-hand side of eq. (2.25). These corrections have been first derived in refs. [12,15]. Moreover, because the color-octet (1) only the part that does contain derivatives. We denote this part with |n; where the fields E 1 and E 2 are the chromoelectric fields computed at the positions x 1 and x 2 , respectively. The contribution of |n; x 1 , x 2 (1) P -wave to the contact term for the color-octet The matrix elements still containing heavy quark and antiquark fields in eq. (3.8) can be evaluated as follows: and k =n The ket state T a |p; where the color indices have been assigned to the states according to footnote 3. Note that necessarily both |p; x 1 , x 2 (0) and |k; x 1 , x 2 (0) have to transform as octet states. Plugging eqs. (3.9) and (3.10) into eq. (3.8) and retaining only terms relevant for P -wave matrix elements, we obtain where the derivatives ∇ 1 and ∇ 2 will eventually act on the wavefunctions once eq. (3.11) is inserted into eq. (2.23). We suppress here and in the following the quark and antiquark positions in the eigenstates when all the positions are the same: The gluonic matrix elements of eq. (3.11) can be cast in the form of an integral over a temporal correlator of chromoelectric fields. We proceed as follows. First, we write k =n In the first equality, we have replaced the ratio of the time independent matrix element to the square of the energy difference in left-hand side with the time integral over the matrix element of the time dependent chromoelectric field supplemented, in a generic gauge, with the Schwinger line in the right-hand side. 5 Note that, whereas the Schwinger lines are path ordered in the color matrices, the fields in the matrix element are time ordered. In the second line of eq. (3.12), the Schwinger lines are in the fundamental representation, whereas, in the third one, the Schwinger line has to be understood in the adjoint representation. 6 The condition k = n in the first and second line of eq. (3.12) has been lifted for r = 0 in the third line, since, according to footnote 4, n ∈ S implies that (0) n|Φ bc (0, x 1 )T c |Ω vanishes after taking the color trace. Because of this, the sum over the states |k (0) is complete and has been replaced with the identity in the last equality of eq. (3.12). Finally, in the last equality we have also computed the color trace: , again using that (0) n| behaves like a color singlet for r = 0. Inserting eq. (3.12) and its Hermitian conjugate in eq. (3.11), we arrive at where the tensor E ij is defined by Since gE e,j 1 (t ) Φ ec 0 (0, x 1 ; t , x 1 ) Φ bc (0, x 1 ) does not contain color matrices, the matrix el- Hence, we could extend the sum in the first line of eq. (3.15) to all n and use the completeness of the eigenstates |n (0) to replace their sum with the identity operator in the second equality: n |n; i, i (0) (0) n; j, j| = δ ij δ ij = N c . Because of translational invariance, we have dropped in the last line of eq. (3.15) the space coordinate, which is the same, but arbitrary, for all operators there. 5 While eq. (3.12) is valid in any gauge, it can be easily verified in the temporal gauge (A0 = 0) where the Schwinger line becomes unity. 6 Consider an infinitesimal time interval dt and fields located at a same but unspecified position, then it holds See also ref. [19]. In eq. (3.15), the chromoelectric field at time t is connected to the origin 0 by the Schwinger line Φ ec 0 (0; t ), which then continues to infinity in the direction. Analogously, the chromoelectric field at time t is connected to the origin 0 by the Schwinger line Φ †ad 0 (0; t), which then continues to infinity in the direction. For a suitable choice of the sign of 0 , the fields in gE e,j (t )Φ ec 0 (0; t )Φ bc are time ordered (T ), and those in Φ †ab Φ †ad 0 (0; t)gE d,i (t) are anti-time ordered (T ). Hence, eq. (3.15) can be interpreted as a cut diagram, which can be useful for perturbative QCD. We show this configuration of Wilson lines graphically in figure 1.
We derive now two properties of the (time ordered) operator gE e,i (t)Φ ec 0 (0; t)Φ bc . First, the Hermitian conjugate of gE e,i (t)Φ ec 0 (0; t)Φ bc is where we have dropped the time ordering prescription whenever the operators appear in the right time or anti-time ordering. This guarantees that the operator Φ †ab Φ †ad 0 (0; t)gE d,i (t) gE e,i (t)Φ ec 0 (0; t)Φ bc is Hermitian. Second, under a gauge transformation U (θ(x)), the operator gE e,i (t )Φ ec 0 (0; t )Φ bc transforms as an adjoint field strength tensor at infinity in the direction: After the first equality we have explicitly required that the fields are time ordered, and in the last equality we have used that U † = U −1 , because U is unitary, and U † = U T , because U is real in the adjoint representation. 7 The gauge transformation property of the operator gE e,i (t )Φ ec 0 (0; t )Φ bc guarantees the gauge invariance of the tensor E ij . Having matched the color-octet contact term, we can compute the color-octet matrix element Ω|O h Q ( 1 S [8] 0 )|Ω from eq. (2.23). Since the product φ (0) 1 P 1 (r)φ (0) * 1 P 1 (r) is isotropic after summing over the polarizations of the 1 P 1 state, the tensor ∇ i r ∇ j r in eq. (3.14) can be replaced by ∇ r · ∇ r δ ij /3. Then, we obtain at leading nonvanishing order in v where E is the dimensionless gluonic correlator The correlator E is the isotropic part of E ij , i.e. (N c δ ij /9) E. The factor 3/N c in the definition of E has been chosen so that eq. (3.18) resembles the pNRQCD expression for the color-octet decay matrix element [14,15] where the correlator E 3 is defined by The contact terms for the operators involving the χ QJ can be computed in a similar way. We obtain where we have again displayed only the terms relevant for the P -wave LDMEs. The tensors T ij 1J are spin projectors that are defined by The wavefunction of a 3 P J state with polarization λ is given by where |1S is a spin-triplet state, and 1M ; 1S|Jλ are Clebsch-Gordan coefficients. At leading order in v, the radial wavefunction R (0) (r) does not depend on the heavy quark spin.
Inserting the above formulas in eq. (2.23), we obtain the following expressions for the LDMEs, which are valid at leading nonvanishing order in v, Together with eqs. (3.5) and (3.18), they enter the NRQCD factorization formulas for production of P -wave quarkonia at leading order in v. Note that Ω|O χ QJ ( 3 S 1 )|Ω depends on E and not on E ij because the matrix element is polarization summed. The results for the color-singlet and color-octet matrix elements lead to the relations These relations are valid for any P -wave quarkonium state, since the right-hand side is independent of the flavor and of the principal quantum number. 8 It follows that all NRQCD matrix elements for production of P -wave quarkonia at leading order in v are fixed once the quantity E and the wavefunctions at the origin are known. By using eqs. (3.5), (3.18), (3.26a), and (3.26b), we obtain the following pNRQCD expressions for the production rates of P -wave quarkonia which are valid at leading order in v and up to corrections of order 1/N 2 c .

Consistency with NRQCD factorization
As we have discussed earlier, the validity of NRQCD factorization requires the shortdistance coefficients to be free of infrared divergences. The short-distance coefficients are determined by replacing the heavy quarkonium state Q in the NRQCD factorization formula (2.1) by a perturbative QQ state, and computing both sides as perturbation series in α s . Hence, for the NRQCD factorization to be valid, the infrared divergences on the lefthand side of eq. (2.1) must be reproduced exactly by the LDMEs on the right-hand side. For our expressions of the LDMEs and perturbative QCD calculations of the short-distance coefficients to be consistent, the pNRQCD expressions of the LDMEs, when computed in perturbative QCD, must reproduce the infrared divergences that appear in the LDMEs J )|Ω scales like 1/m 2 times a constant that is independent of the flavor and the principal quantum number has been assumed on phenomenological grounds in ref. [31].
when the quarkonium states are replaced by perturbative QQ states. This is nontrivial in the case of the color-octet LDMEs, which involve arbitrary exchanges of gluons between the quark, antiquark, and the gauge-completion Wilson lines.
Well-established methods to investigate infrared divergences in perturbative QCD calculations have been developed in proofs of perturbative factorization [32][33][34][35][36]. This consists of applying approximations to gluon propagators and vertices that simplify the perturbative QCD expressions while preserving the infrared divergences, and reorganizing the resulting expressions in terms of Wilson lines. In ref. [19], the authors considered the infrared divergences in the color-octet LDME Ω|O Q(p 1 )Q(p 2 ) ( 3 S [8] 1 )|Ω . In figure 2, we show some representative Feynman diagrams that appear in the perturbative QCD calculation of the LDME. The quark and antiquark in the final state are on shell and have momenta p 1 = p+q and p 2 = p − q, respectively. There are additional diagrams that can be obtained by moving gluon attachments from one side of the cut to the other side, or by replacing gluon attachments to the quark (antiquark) line by attachments to the antiquark (quark) line. Note that, since the scale m has been integrated out in NRQCD, gluons cannot have momenta that exceed scales of order mv. In these diagrams, infrared divergences arise when the virtual gluons have soft momenta, that is, when the temporal and spatial components of the gluon momenta are small and of the same order. In such case, the soft approximation can be used for propagators and vertices involving these gluons [33,37,38]. When the quark and antiquark in the final state are on shell, the soft approximation leads to the eikonal approximation, where a soft gluon with momentum k has propaga-tor and vertex given in Feynman gauge by i/(β · k + iε) and ∓ig s T a β µ , where β is the four-velocity of the quark or antiquark line, and the minus (plus) sign applies to the quark (antiquark) line. Note that the vertex factor in the eikonal approximation does not involve gamma matrices. It is clear that the Feynman rules for eikonal gluon attachments to the quark line are equivalent to the Feynman rules for gluon attachments to the path-ordered Wilson line in the fundamental representation in the direction of the quark momentum. Similarly, eikonal gluon attachments to the antiquark line are equivalent to gluon attachments to anti-path-ordered Wilson line in the direction of the antiquark momentum. That is, the eikonal approximation decouples soft gluon attachments to quark and antiquark lines, and factors them out in the form of Wilson lines. If the QQ in the final state is in a color-singlet state, we obtain the following infrared factor by collecting all Wilson-line factors arising from the application of the eikonal approximation to soft gluon attachments to the quark and antiquark lines on both sides of the cut: where Φ p 1 = Pexp −ig s 1 )|Ω /I(p 1 , p 2 ) is free of such divergences. This infrared factor has first been obtained in Refs. [18,19] in the analysis of gluon fragmentation functions for the process g → QQ + X.
The quantity Ω|O Q(p 1 )Q(p 2 ) ( 3 S 1 )|Ω /I(p 1 , p 2 ) is not completely free of infrared divergences, because the LDME can have additional infrared divergences from gluon exchanges between the quark and antiquark lines. These infrared divergences are not captured by the infrared factor (3.29), because they arise from the region of gluon momentum where its temporal component is much smaller than the spatial components; that is, they are singularities associated with potential gluons. From proofs of perturbative factorization in quarkonium decays in NRQCD, it can be seen that at leading order in v infrared divergences arising from gluon exchanges between the quark and antiquark in a color-singlet state can be absorbed into the quarkonium wavefunctions at the origin [5]. Logarithmic divergences arising from quarkonium wavefunctions at the origin appear from two loops, and hence, these divergences affect the color-octet LDME from three loops, because at least one extra gluon is needed in order to carry the color charge over the cut when the QQ is in a color-singlet state. The Feynman diagram in figure 2(f) shows a contribution at NNNLO in α s from which such logarithmic divergences appear. It is possible that there are additional infrared divergences that cannot be absorbed into the quarkonium wavefunctions at the origin nor the infrared factor, if there are infrared divergences arising from soft gluons interacting with potential gluons. Nonetheless, such effects do not appear at least until three-loop accuracy, as can be seen in explicit calculations of the color-octet LDME in perturbative QCD [20,21]. Therefore, the infrared factor does contain all infrared divergences appearing in the color-octet LDME to two-loop accuracy, which come from soft gluons carrying momenta of order mv.
It is particularly interesting to investigate the infrared factor (3.29) at the lowest nonvanishing power in the relative momentum q, because matching calculations are usually carried out at a specific accuracy in q/m. The lowest power term is of quadratic order in q, where a factor of q comes from each side of the cut. To relative order q 2 /m 2 , the infrared factor is given by I 2 (p, q), which is defined in eq. (32) of ref. [19] as 9 is an adjoint Wilson line along p. Since in eq. (3.30) a momentum q comes from each side of the cut in the squared amplitude, the infrared factor contributes to the production of a color-singlet P -wave state. The same infrared factor appears in the fragmentation of a gluon into color-singlet QQ. Since the gluon fragmentation process produces a QQ in a color-octet 3 S 1 state, the divergences must match the infrared divergences in the color-octet LDME 1 )|Ω , when Q is replaced by a color-singlet QQ state. This agreement has been confirmed explicitly through one-loop and partial two-loop calculations of the coloroctet LDME in ref. [20]; the two-loop calculations have only been done for the diagrams involving the gauge-completion Wilson lines. Since the LDME Ω|O Q ( 3 S [8] 1 )|Ω appears in the NRQCD factorization formula at leading order in v, the same infrared divergences have to occur in the operator Ω|O Q ( 1 S [8] 0 )|Ω because of the heavy quark spin symmetry. In the rest frame of the QQ, where p = 0 and q 0 = 0, Φ p (λ) is the Schwinger line Φ 0 (0, 0; t, 0) in the SU(3) adjoint representation with t = p 2 λ, and p µ q ν G a νµ (λp) is − p 2 q i E a i (t) with the chromoelectric field located at 0. Therefore, the infrared factor I 2 (p, q) defined in eq. (3.30) can be written in the rest frame of the QQ as  (3.14) and (3.23), respectively, when written in momentum space for a color-singlet QQ state with relative momentum q. The pNRQCD expressions for the color-octet LDMEs in eqs. (3.18) and (3.26b) are, therefore, expected to reproduce the same infrared divergences obtained from the NRQCD factorization. This is straightforward to check at one-loop accuracy. By computing the correlator E at order α s accuracy in dimensional regularization with d = 4 − 2 spacetime dimensions, we obtain is the Casimir of the fundamental representation of SU(3) and the subscripts UV and IR indicate the origin (ultraviolet and infrared, respectively) of the 1/ poles. After renormalizing the UV divergence at the scale Λ, E satisfies the following renormalization group equation which, in turn, implies the following renormalization group equation for the NRQCD matrix elements (see eqs. (3.26a) and (3.26b)) The same evolution equation 1 ) . Equation (3.35) agrees with the evolution equation derived from a perturbative calculation in NRQCD [5], and, therefore, the UV divergence at one-loop accuracy of the color-octet LDME in pNRQCD is the same as in NRQCD. Since loop corrections to NRQCD matrix elements are scaleless, UV poles cancel IR poles in eq. (3.33). Hence, the IR divergence at one-loop accuracy of the color-octet LDME in pNRQCD is the same as in NRQCD too.
At two loops, some consistency checks between the pNRQCD result and the NRQCD factorization can be made based on the two-loop calculations in refs. [19] and [20]. In ref. [19], two-loop corrections to the infrared factor I 2 (p, q) associated with the gaugecompletion Wilson lines were computed. These contribute to the infrared divergences of the LDMEs O h Q ( 1 S [8] 0 ) and O χ QJ ( 3 S [8] 1 ) at order α 2 s , and were reproduced in ref. [20] through the explicit calculation of the LDMEs. Since, as we have seen, the calculation of the infrared factor I 2 (p, q) is equivalent to the calculation of the infrared divergences in the contact terms V O( 1 S [8] 0 ) and V O( 3 S [8] 1 ) , and eventually in the chromoelectric field tensor E ij , it follows that the pNRQCD expressions for the color-octet LDMEs have the same infrared divergences associated with the gauge-completion Wilson lines as those found in the NRQCD calculations of ref. [20].
It may be interesting to note that at order α s the pole structure of the correlator E and of the correlator E 3 , defined in eq. (3.21) and relevant for the decay widths of P -wave quarkonia, is the same, i.e. the one given in eq. (3.33). This reflects at the pNRQCD level the fact that at the NRQCD level the one-loop renormalization group equation (3.35) is the same as the one-loop renormalization group equation for the decay matrix elements appearing in inclusive decays of P -wave quarkonia [5]. The observation, however, ceases to hold at two loops, because at this order, E receives contributions from the gauge-completion Wilson lines, which are absent in E 3 . 10 An important issue in NRQCD factorization is whether the color-octet LDMEs are independent of the direction of the gauge-completion Wilson lines, which is necessary in establishing the universality of the NRQCD production matrix elements. While a general argument for the universality has been suggested in ref. [39], an explicit verification has only been done at two-loop accuracy [19][20][21][22]. In the pNRQCD expression of the color-octet LDMEs, the dependence on the direction of the gauge-completion Wilson lines is encoded in the tensor E ij . For the case of polarization-summed cross sections, where the polarization of the quarkonium in the final state is summed over, only the isotropic part of E ij , given by (N c δ ij /9) E, contributes to the color-octet LDMEs, and, therefore, the dependence on the direction of the gauge-completion Wilson lines disappears due to rotational symmetry. Hence, the pNRQCD expressions of the color-octet LDMEs support the universality of the NRQCD LDMEs for polarization-summed cross sections of P -wave quarkonia.
4 Phenomenology of inclusive production of P -wave quarkonia

Phenomenological determination of E
In order to compute the cross sections of P -wave quarkonia, it is necessary to obtain values for the radial wavefunctions at the origin and the gluonic correlator E. Radial wavefunctions at the origin can be computed by solving the Schrödinger equation or can be determined from decay rates. The gluonic correlator E could, in principle, be computed from lattice QCD. Since, however, a lattice QCD determination of E is not available at the time of this study, we determine E phenomenologically by comparing the measured ratio of Pwave quarkonium cross sections to the expression obtained from the NRQCD/pNRQCD factorization.
We use the ratio r 21 of the χ c2 (1P ) and χ c1 (1P ) differential cross sections defined by where p T is the transverse momentum of the χ cJ (1P ). This ratio has been measured at the LHC by CMS [40] and ATLAS [41]. In order to compute the cross sections σ χ cJ (1P ) in NRQCD, we employ the short-distance coefficients σ QQ( 3 P [1] J ) and σ QQ( 3 S [8] 1 ) that were computed in ref. [42] at next-to-leading order (NLO) in α s for a center-of-mass energy of 7 TeV with rapidity range |y| < 0.75. The calculation in ref. [42] used for the charm quark mass m c = 1.5 GeV, CTEQ6M parton distribution functions at the scale µ F = p 2 T + 4m 2 c , and computed α s at the same scale µ R = p 2 T + 4m 2 c , running at two loops with n f = 5 light quark flavors and Λ J ) depend on the scheme and scale Λ used to renormalize the color-octet matrix 10 The difference between E and E3 signals also the violation of the crossing symmetry beyond leading order in the corresponding NRQCD matrix elements [5].
1 )|Ω , which in pNRQCD we identify with the renormalization scheme and scale used for E. We renormalize in the MS scheme at the scale Λ = m c , where, as above, m c = 1.5 GeV is the charm quark mass. 11 We estimate the uncertainties in the short-distance coefficients to be 30% of their central values, which account for corrections of relative order α s or v 2 that we neglect. The variation of the scale µ F for the parton distribution functions and of the renormalization scale µ R for α s affects the short-distance coefficients by less than 25% of the central values.
We neglect the uncertainty of order 1/N 2 c in the wavefunctions (see eq. (2.21)) compared to other uncertainties. We use the pNRQCD expressions for the matrix elements in eqs. (3.26). The quarkonium wavefunctions at the origin cancel in the ratio r 21 , which makes r 21 an attractive observable to extract E.
Because the transverse momentum p T of the quarkonium can be much larger than m c in the kinematical range that we consider, the logarithms of p T /m c can have a large impact on the perturbative corrections to the short-distance coefficients. The resummation of the logarithms of p T /m c have been computed in refs. [42,43] at leading logarithmic accuracy at leading power (LP) in the expansion in powers of m c /p T . Following refs. [42,43], we use the label LP+NLO for the short-distance coefficients including the additional LP contributions that augment the fixed-order NLO calculations. The results in ref. [42,43] show that the additional LP contributions are numerically significant. However, we note also that the large additional LP corrections come mainly from partial contributions of order α 2 s , and may overestimate the size of the corrections to the short-distance coefficients at nextto-next-to-leading order in α s . Therefore, in computing P -wave charmonium production cross sections, we consider both the fixed-order NLO and LP+NLO calculations of the short-distance coefficients, which can provide an estimate of the uncertainty coming from uncalculated corrections of higher orders in α s .
In order to compare to measurements, we compute the values of r 21 multiplied with B χ c2 (1P ) /B χ c1 (1P ) , where B χ cJ (1P ) = Br(χ cJ(1P ) → J/ψγ) × Br(J/ψ → µ + µ − ), and Br stands for the branching ratio. We compute B χ cJ (1P ) from the PDG values [44]. Since the measurements of r 21 are given as functions of the transverse momentum p J/ψ T of the J/ψ, we compute p J/ψ T from the transverse momentum p T of the χ cJ (1P ) using which is a good approximation since m J/ψ ≈ m χ cJ (1P ) . Corrections to eq. (4.2) affect the cross section by less than 1% in the kinematical range that we consider [42]. By performing a least-squares fit to the measured values of r 21 × B χ c2 (1P ) /B χ c1 (1P ) by CMS [40] and ATLAS [41], we obtain, from fixed-order NLO calculations of the short-distance coefficients, E| NLO (Λ = 1.5 GeV) = 1.17 ± 0.05, (4.3) 11 Although from an effective field theory point of view it would make sense to choose for Λ a scale close to the soft scale mv, the specific choice of scale for Λ is without consequences here, since, according to the renormalization group equation (3.34), the Λ dependence of E is encoded into an additive constant that may be freely reshuffled between E and the short distance coefficients. 3) and (4.4) reflects the difference between the fixed-order NLO and LP+NLO calculations of the shortdistance coefficients. We show our result for r 21 compared to ATLAS and CMS data in figure 3. In the following sections, we will use both values of E in eqs. (4.3) and (4.4) to compute cross sections of χ cJ and χ bJ states at the LHC. More precisely, when computing cross sections of χ cJ , we will use the result in eq. (4.3) with the fixed-order NLO expressions of the short-distance coefficients, and we will use the result in eq. (4.4) with the LP+NLO expressions of the short-distance coefficients. When computing cross sections of χ bJ , we will combine the two determinations of E into E(Λ = 1.5 GeV) = 2.8 ± 1.7,

Production of χ cJ (1P )
We compute the inclusive production cross sections of χ cJ (1P ) from proton-proton collisions at the LHC based on the expression of the matrix elements given in eqs. (3.26). We use the same fixed-order NLO and LP+NLO short-distance coefficients from ref. [42] that we used in section 4.1, and use the corresponding determinations of E at the scale Λ = 1.5 GeV given in eqs. (4.3) and (4.4) for fixed-order NLO and LP+NLO, respectively. We determine the value of the charmonium 1P -state wavefunction at the origin from the two-photon decay rates of the χ c0 (1P ) and χ c2 (1P ). For consistency with our calculation of the cross sections, we use the NRQCD factorization formulas for the decay rates at leading order in v, while we include order α s corrections to the short-distance coefficients at the amplitude level. The pNRQCD expressions for the two-photon widths at leading order in v read [5,15,17] where e c = 2/3 is the fractional electric charge of the charm quark, and α is the QED coupling constant. In the decay rates, we use α = 1/137 reflecting the fact that the photons in the final states are on shell, and use α s = 0.282, which is evaluated at the scale m χ cJ /2. We use m c = 1.5 GeV for consistency with the calculation of the cross sections. By comparing these formulas with the BESIII measurements of the decay rates [45], we obtain |R (0) (0)| 2 = 0.041 GeV 5 for the χ c0 (1P ), and |R (0) (0)| 2 = 0.073 GeV 5 for the χ c2 (1P ).
For the calculation of the NRQCD matrix elements, we take the average |R (0) (0)| 2 = 0.057 GeV 5 as the central value. We attribute to the central value a 30% uncertainty, which accounts for the uncalculated corrections of relative order v 2 . We note that, while the order-α 2 s corrections to the two-photon widths have been computed in ref. [46], these corrections depend strongly on the scheme and the scale associated with the renormalization of the wavefunctions at the origin (see, for example, refs. [47,48]). The strong dependence on the scheme and scale will not cancel in production rates unless the short-distance coefficients for cross sections are also computed to the same accuracy in α s . Since the short-distance coefficients for heavy quarkonium production are currently known at NLO accuracy in α s , we neglect the order-α 2 s corrections to the two-photon widths. The results for the χ c1 (1P ) and χ c2 (1P ) cross sections are shown in figure 4 against ATLAS data. The fixed-order NLO and LP+NLO calculations of the short-distance coefficients give results for the cross sections that are compatible within uncertainties. We note that the NLO results for the cross sections are systematically lower than the ATLAS data in particular at low p T , while the LP+NLO results agree well with the measurements.

Polarization of χ cJ (1P )
In the case of polarized cross sections, the non-isotropic part of E ij can in principle contribute to the color-octet matrix elements, and, if such contribution is nonvanishing, the color-octet matrix elements can acquire a dependence on the arbitrary direction of the gauge-completion Wilson lines. For the universality of the NRQCD matrix elements to be valid also for the case of polarized cross sections, such non-isotropic contributions should vanish in the NRQCD matrix elements. This has not been proved. However, it is common practice in the phenomenological literature to assume that the polarized color-octet LDMEs do not depend on the direction of the gauge-completion Wilson lines, and, therefore, do not depend on the polarization λ (see for instance refs. [42,[49][50][51][52][53]). This amounts at assuming that the polarized color-octet LDMEs Ω|χ † σ i T a ψΦ †ab (0)P χ QJ (λ,P =0) Φ bc (0) ψ † σ i T c χ|Ω can be related for any polarization λ of the χ QJ state to the polarizationsummed color-octet LDME Ω|O χ QJ ( 3 S [8] 1 )|Ω , which is independent of the direction of the gauge-completion Wilson lines, through the relation (4.8) and similarly for h Q : (4.9) These relations are consequences of the heavy quark spin symmetry and the assumption of the universality of the NRQCD matrix elements for polarized quarkonia in the right-hand sides of eqs. (4.8) and (4.9). Similar relations hold for the color-singlet LDMEs: In this case, however, they are not assumptions, but follow from the vacuum-saturation approximation and rotational symmetry. Under the universality assumption of the polarized color-octet LDMEs, we can compute the polarization of χ c1 (1P ) and χ c2 (1P ) produced at the LHC. The polarization parameters λ χ c1 θ and λ χ c2 θ are defined by where ξ χ cJ is the fraction of J/ψ produced with longitudinal polarization from decays of χ cJ (1P ). We use the hadron helicity frame to define the spin quantization axis of the J/ψ. The polarized cross sections can be computed by using the short-distance coefficients from ref. [42] for the polarized production of χ cJ (1P ), and eq. (4.8) to relate the color-octet LDMEs for polarized and polarization summed χ cJ (1P ) states. Equation (4.8) relies on the assumption that the NRQCD matrix elements are universal for polarized cross sections.
Our results for λ χ c1 θ and λ χ c2 θ are shown in figure 5 against experimental constraints coming from the CMS data [54].
It has been pointed out in refs. [53,[55][56][57] that modifications to the central values of χ cJ cross section measurements in refs. [40,41] due to polarizations of χ c1 and χ c2 may be important. The central values in refs. [40,41] have been obtained by assuming isotropic decay angular distributions, while the scale factors that modify the central cross-section values are listed for extreme polarization scenarios. We compute the scale factors by using our results for λ χ c1 θ and λ χ c2 θ in figure 5, and by assuming that the scale factor depends smoothly on polarization. The scale factors for r 21 that we obtain are about 0.95 for CMS data [40], and about 0.89 for ATLAS data [41]. These scale factors have mild effects on our phenomenological determinations of E, and slightly reduce the central values by less than the uncertainties in E. In ref. [53], a smaller scale factor of 0.85 has been obtained by using a data-driven analysis. In this case, the values of E determined from r 21 reduce by slightly more than the estimated uncertainties. Nonetheless, the changes in the central values of E have negligible effect on our phenomenological results for P -wave charmonia and bottomonia.

Production of χ bJ (nP )
Owing to the universality of the correlator E, we can compute in pNRQCD, at leading order in v, inclusive production cross sections of χ bJ (nP ) from proton-proton collisions at the LHC without having to fit new octet LDMEs. We first consider the ratio r 21 for the differential cross sections of χ b2 (nP ) and χ b1 (nP ) states. We compute the short-distance coefficients at next-to-leading order accuracy using the FDCHQHP package [58]. We use the bottom quark mass m b = 4.75 GeV, CTEQ6M parton distribution functions at the scale µ F = p 2 T + 4m 2 b , and compute α s at the same scale running at two loops with n f = 5 light quark flavors and Λ (5) QCD = 226 MeV. Since the range of p T for the χ bJ (nP ) that we consider is not too large compared to the mass of the χ bJ (nP ) states, we do not resum logarithms of p T /m b . We take the renormalization scale Λ of the color-octet matrix element Ω|O χ bJ (nP ) ( 3 S [8] 1 )|Ω , which in pNRQCD is the renormalization scale of E, to be m b . We compute E at the scale 4.75 GeV by using the one-loop renormalization-group-improved formula where β 0 = 11N c /3 − 2n f /3 and E(m c ) is given in eq. (4.5). We take into account the uncertainty in E, and estimate the uncertainties from uncalculated corrections of order v 2 to be 10% of the central values.
Our result for r 21 in the case of χ bJ (1P ) states is shown in figure 6 against LHCb [59] and CMS [60] data. The result is given as a function of the transverse momentum p In figure 6, we did not include the feeddown contributions, which come dominantly from decays of Υ(2S) and Υ(3S). Since the branching ratios into χ b1 (1P ) and χ b2 (1P ) are almost the same, and our results for r 21 is close to unity, we expect our results for r 21 to be almost unchanged by the inclusion of the feeddown contributions. In pNRQCD, the wavefunction dependence factorizes in the cross sections, see eqs. (3.28), and cancels at leading order in v when considering the ratio of states with the same principal quantum number. Moreover, the ratio r 21 for the χ bJ (1P ) states in figure 6 is almost independent of p T . It is, therefore, a specific prediction of (leading order) pNRQCD that the ratio r 21 is independent of the principal quantum number, i.e. that the ratios for the χ bJ (2P ) and χ bJ (3P ) states is expected to be the same as the one shown in figure 6 for the χ bJ (1P ) states.
For the computation of the production rates of χ bJ (nP ) we need the values of the P -wave bottomonium wavefunctions at the origin. Since the decay widths of the χ bJ (nP ) states are in general not well known, differently from the χ cJ (1P ) states, we take the central values of the wavefunctions at the origin to be the average of the potential-model calculations considered in ref. [17]. They are   Figure 7. Differential production cross sections of the χ b1 (nP ) and χ b2 (nP ) (n = 1, 2, and 3) at the LHC center of mass energy √ s = 7 TeV and in the rapidity range 2 < y < 4.5. The χ b1 (nP ) cross sections have been multiplied by 10 to distinguish them from the χ b2 (nP ) cross sections on the logarithmic scale.
While the p T -dependent absolute production rates of χ bJ states at the LHC have not been reported yet, we compute cross sections of χ b1 (nP ) and χ b2 (nP ) (n = 1, 2, and 3) at the LHC with √ s = 7 TeV in the rapidity range 2 < y < 4.5. Our results for the differential cross sections are shown in figure 7. Once measurements become available, our predictions could be compared with data.
To compare with data we need to compute the feeddown fraction R χ b (nP ) Υ(n S) , which is defined as the fraction of Υ(n S) produced from decays of χ b1 (nP ) and χ b2 (nP ). We compute R χ b (nP ) Υ(n S) by using the formula where the inclusive cross sections σ χ bJ (nP ) and σ Υ(n S) are the sum of the direct production rate and the feeddown contribution. We neglect the contributions from χ b0 (nP ), because the branching ratios Br(χ b0 (nP ) → Υ(n S) + γ) are small, while the cross sections σ χ b0 (nP ) are similar in size compared to σ χ b1 (nP ) or σ χ b2 (nP ) . The feeddown contribution to the Υ(n S) production rate is given by where we truncate the sum at n = 3. In the production rate of Υ(3S), we only consider the feeddown contribution from decays of χ b (3P ). In the production rate of Υ(2S), we consider the feeddown contributions from decays of χ b (3P ) and χ b (2P ), and also from decays of Υ(3S). For the Υ(1S) production, we consider the feeddowns from decays of , Υ(3S) and Υ(2S) states. For the production of the χ bJ (3P ) state we neglect the contribution from feeddowns, while the feeddown contribution for the χ bJ (nP ) production rate for n = 1, 2 is given by where we truncate the sum again at n = 3. While we compute the direct production cross sections of χ bJ (nP ) in pNRQCD from eq. (3.28b) (the result has been shown in figure 7), we compute the direct production rates of Υ(n S) from the NRQCD factorization formula: which is truncated at relative order v 4 . The NRQCD operators for production of Υ(nS) are defined by (4.19d) We have used the heavy quark spin symmetry to relate the LDMEs Ω|O Υ(n S) ( 3 P J ) by inserting color matrices and gauge-completion Wilson lines between the quark and antiquark fields. As is commonly done in phenomenological studies of S-wave heavy quarkonium production, we neglect contributions from color-singlet matrix elements of relative orders v 2 and higher that are obtained by inserting powers of ← → D 2 between the quark and antiquark fields, because their contributions to the cross section turn out to be numerically small. Since the coloroctet matrix elements for the Υ(n S) states have not been computed from first principles, we take the values that were determined in ref. [61] from fits to measured cross sections at the LHC. We take the measured values of the branching ratios from ref. [44]. Since the branching ratios of the χ bJ (3P ) states have not been measured yet, we use the theoretical prediction of ref. [61]. For the branching ratios Br(Υ(3S) → Υ(1S) + X) and Br(Υ(2S) → Υ(1S) + X), we take the sum of the branching ratios into Υ(1S)π + π − and Υ(1S)π 0 π 0 . Our pNRQCD results for the feeddown fractions R χ b (nP ) Υ(n S) are shown in figure 8 against LHCb data [62]. Compared to the previous pNRQCD results in ref. [23], the results here have larger uncertainties because we consider a wider range of values for E, which reflect the uncertainty due to unknown corrections of higher orders in α s in the short-distance coefficients. Also, the results here include feeddown contributions to χ b (2P ) and χ b (1P ) cross sections, which increase slightly the fractions R χ b (2P ) Υ(n S) and R χ b (1P ) Υ(n S) , and improve the agreement with the LHCb data.

Conclusions
In this work, we have developed a formalism for computing the NRQCD long-distance matrix elements for the inclusive production of heavy quarkonia, based on pNRQCD. Our formalism applies to strongly coupled quarkonia, which include excited charmonium and bottomonium states. We obtain expressions of the LDMEs that are given by products of quarkonium wavefunctions at the origin and universal gluonic correlators, analogously to the pNRQCD expressions of decay LDMEs. The computation of the production LDMEs in this paper is valid up to corrections of relative order 1/N 2 c and v 2 . Corrections of higher orders in v come from higher order corrections in the pNRQCD expressions of the NRQCD LDMEs, as well as from corrections to the quarkonium wavefunctions.
Based on the general formalism that we have developed in section 2, we have computed NRQCD LDMEs for hadroproduction of P -wave heavy quarkonia at leading nonvanishing order in v in section 3. For the case of color-singlet LDMEs, we reproduce the known result from the vacuum-saturation approximation where the production LDMEs and decay LDMEs are both given by the modulus square of the derivative of the quarkonium wavefunctions at the origin. For the color-octet LDMEs, we obtain expressions that are similar to the color-octet decay LDMEs, except that they depend on a gluonic correlator E that is different from the one found for decay LDMEs, due to a different field arrangement and the gauge-completion Wilson lines that are necessary to ensure the gauge invariance of the color-octet LDMEs. The results confirm our previous calculation of the P -wave LDMEs in ref. [23], where we assumed heavy-quark spin symmetry to reduce the number of independent LDMEs. In this work, we explicitly confirm the validity of the heavy-quark spin symmetry for P -wave LDMEs at leading nonvanishing order in v.
Our results for the color-octet P -wave LDMEs help identify the infrared behavior, which is necessary in testing the validity of the NRQCD factorization in heavy quarkonium production. Since the color-octet P -wave LDMEs are given by products of the modulus square of the derivatives of the wavefunctions at the origin and the gluonic correlator E, defined in eq. (3.19), they have additional infrared divergences that come from E compared to the color-singlet P -wave LDMEs. These are consistent with the infrared factor found computing the gluon fragmentation function in ref. [19], which, in turn, agrees with the explicit calculation of the infrared divergences in the color-octet LDMEs in ref. [20]. Therefore, our results are consistent with the known infrared behavior of the NRQCD LDMEs. Furthermore, our expressions for the color-octet P -wave LDMEs for polarization-summed cross sections are independent of the direction of the gauge-completion Wilson lines. This implies that the color-octet P -wave LDMEs are process independent, which is a necessary condition for the validity of the NRQCD factorization. Hence, our results for the P -wave LDMEs support the validity of the NRQCD factorization in heavy quarkonium production for polarization-summed cross sections.
We have computed production cross sections of χ cJ and χ bJ states at the LHC for J = 1, 2 based on our calculations of the P -wave LDMEs and the NRQCD short-distance coefficients that are available through next-to-leading order accuracy in α s . Since a lattice QCD determination of the correlator E is not available, we determine E by comparing the theoretical expression of the differential cross section ratio (dσ χ c2 (1P ) /dp T )/(dσ χ c1 (1P ) /dp T ) with measurements at the LHC by CMS [40] and ATLAS [41]. This allows a phenomenological determination of E that does not depend on the value of the wavefunction at the origin. This determination improves our previous determination in ref. [23] where absolute cross section measurements were used. After fixing the value of the derivative of the P -wave charmonium wavefunction at the origin on the measured two-photon decay rates, we have computed the cross sections of χ c1 (1P ) and χ c2 (1P ) at the LHC, which agree well with ATLAS measurements [41]. Our results for the χ c1 (1P ) and χ c2 (1P ) cross sections are consistent with our previous results in ref. [23]. We have also computed the polarizations of χ c1 (1P ) and χ c2 (1P ) at the LHC under the assumption that the NRQCD factorization holds for polarized cross sections. Our results are consistent with experimental constraints from CMS [54]. The universality of the correlator E allows to compute P -wave bottomonium cross sections without having to fit any new octet LDME. From a phenomenological point of view, this is the most relevant gain in the pNRQCD approach. In the bottomonium case, we have computed the differential cross section ratio (dσ χ b2 (1P ) /dp T )/(dσ χ b1 (1P ) /dp T ) and the feeddown fractions R χ b (nP ) Υ(n S) at the LHC, for which measurements are available. We find good agreements with data [59,60,62]. We have improved our results for the feeddown fractions R χ b (nP ) Υ(n S) compared to our previous results in ref. [23] by considering also the contributions from cascade decays. This results in a better agreement with data. Finally, we made predictions for absolute production rates of χ b1 (nP ) and χ b2 (nP ) (n = 1, 2, and 3) at the LHC, which may be compared with data once measurements become available.
The phenomenological results in this work are based on the determination of the correlator E. We have obtained it from measured cross sections of charmonia. It would be desirable, however, to have a lattice QCD determination of this nonperturbative quantity. Apart from its phenomenological relevance, a lattice QCD determination of E would also provide a nontrivial test of QCD in the nonperturbative regime. While in the charmonium sector, the derivative of the wavefunctions at the origin can be extracted from two-photon decay data, this is not yet possible in the bottomonium sector for lack of data. In the meantime, relying on models is a major source of uncertainty, and reducing this uncertainty a crucial step to further progress. Some advancements towards an accurate determination of the quarkonium wavefunctions at the origin also for strongly coupled quarkonia has been made recently in refs. [47,48]. From a general perspective, it is important to emphasize that our phenomenological results, having been derived in an effective field theory framework, can be systematically improved through the formalism developed in this work by including corrections of higher order in v. These come from higher dimensional operators in the NRQCD factorization formula, from higher order corrections to the pNRQCD expansion of the NRQCD long-distance matrix elements, and from higher order corrections to the wavefunctions originating from higher order corrections to the pNRQCD potential.
Higher order corrections in α s can be included within perturbative QCD. Studies of relativistic corrections to the NRQCD factorization formula from higher dimensional operators in refs. [63,64] suggest that corrections of higher orders in v can be appreciable for inclusive production of J/ψ, although phenomenological applications have been limited so far by the lack of determinations of these higher order long-distance matrix elements.
As a first and obvious outlook, we expect the formalism developed in this work for computing production LDMEs to be applicable to any strongly coupled quarkonia, not only the P -wave quarkonia considered here and in [23]. These include the important case of S-wave quarkonia. Calculations of the S-wave quarkonium production LDMEs may shed light on the longstanding puzzle of the J/ψ production mechanism (as long as the charmonium ground state may be assimilated to strongly coupled quarkonia), as well as on issues in the polarization of J/ψ and Υ(nS), and the η c production rate. Further, it would be interesting to extend the formalism to describe the production of weakly coupled quarkonia, which are possibly the lowest-lying quarkonium states, and also the production of quarkonium exotica (hybrids, tetraquarks). In these cases, the formalism will require the introduction of nonperturbative matrix elements or functions that will be different from the ones encountered in the strong coupling regime. Finally, we also expect this approach to be useful in studying heavy quarkonium production in heavy ion collisions, and help in this way to unveil the nature of the hot and dense phase of QCD.