Resummation and renormalization of kinematical effects in inclusive $P$-wave quarkonium production

We investigate the renormalization properties of the shape function formalism for inclusive production of $P$-wave heavy quarkonia, which arises from resumming a class of corrections coming from kinematical effects associated with the motion of the heavy quark and antiquark pair relative to the quarkonium. Such kinematical effects are encoded in the nonperturbative shape functions, which are normalized to the corresponding nonrelativistic QCD long-distance matrix elements. By using the known ultraviolet divergences in the matrix elements, we derive the large-momentum asymptotic behavior of the shape functions. This strongly constrains the form of the shape functions and significantly reduces the dependence on the nonperturbative model. Based on these results we show that the shape function formalism at loop level can be useful in taming the threshold logarithms at large transverse momentum, and at small transverse momentum the kinematical corrections reduce the sizes of $\chi_c$ and $\chi_b$ cross sections which may improve agreement with measurements.


Introduction
First-principles based studies of inclusive heavy quarkonium production phenomenology have usually been done in the nonrelativistic QCD (NRQCD) factorization formalism [1]. In this formalism, the production rate of a heavy quarkonium is factorized into products of perturbatively calculable short-distance coefficients and nonperturbative long-distance matrix elements. The predictive power of the factorization formalism comes from the universality and the nonrelativistic power counting of the matrix elements, which allow quantitative descriptions of production rates once values of a few nonperturbative parameters are determined from experiments. This approach had limited success in deciphering the heavy quarkonium production mechanism: as have been discussed in the literature, results from the global fit analyses of J/ψ production based mostly on total inclusive production and low transverse momentum data [2][3][4][5][6][7] are in conflict with analyses mostly based on large transverse momentum data [8][9][10][11][12][13][14][15][16] and with LHC measurements of polarized and polarizationsummed hadroproduction rates at large transverse momentum. Similarly, in the case of ψ(2S), NRQCD analyses generally have difficulty explaining production rates with transverse momentum smaller than 2 -3 times the heavy quarkonium mass [10,11,14,16,17]. Moreover, some matrix element determinations can lead to large uncertainties in the cross sections at very large transverse momentum due to mixings of contributions from different channels, which comes from renormalization of the long-distance matrix elements [17]. These call for investigation of possible refinements to the NRQCD factorization approach in quarkonium production phenomenology that can modify the transverse momentum dependences of the cross sections.
In an inclusive production process, the quark Q and antiquarkQ produced in a hard process can have nonzero total momentum relative to the heavy quarkonium. This happens because the QQ can emit an arbitrary number of soft particles before evolving into a heavy quarkonium, which carry small but nonzero momenta. This kinematical effect gives rise to operators that are given in the form of total derivatives of the QQ bilinear [18][19][20]. While these operators are suppressed by powers of the heavy-quark velocity v in the nonrelativistic power counting of NRQCD, they can be enhanced due to threshold effects associated with the boundary of phase space [19,20]. The corrections from these operators can be partially resummed to all orders in v, which leads to the shape function formalism developed in ref. [21], where products of short-distance coefficients and NRQCD matrix elements are replaced by convolutions of short-distance coefficients and nonperturbative shape functions. Here, shape functions are functions of the momentum of the QQ pair in the rest frame of the heavy quarkonium in the final state. A similar formalism based on the soft-collinear effective theory have been developed in ref. [22]. It has been anticipated that the kinematical effects computed in the shape function formalism will have a significant impact on the transverse-momentum differential cross sections of heavy quarkonia. However, phenomenological applications of this formalism have been limited, because shape functions are nonperturbative functions that are in general unknown, other than the fact that they are normalized to the corresponding NRQCD matrix elements. That is, unless one can determine the shape functions from first principles, the shape function formalism has less predictive power than the NRQCD factorization formalism.
Another shortcoming of the shape function formalism is that, in its original form, it does not correctly incorporate the renormalization of the NRQCD matrix elements. Because the shape functions are normalized to the NRQCD matrix elements, their normalizations must reproduce the same ultraviolet (UV) divergences in the matrix elements for the formalism to be consistent with NRQCD. This point has not been considered in the development of the formalism in refs. [21,22], because at that time, phenomenological applications of the NRQCD factorization formalism have usually been done at leading order in the strong coupling.
The UV divergences in the normalizations of the shape functions can have important phenomenological implications. First, the known UV divergences fix the large-momentum asymptotic behaviors of the shape functions, which strongly constrain their forms. This happens because, as we will see later, the UV divergence comes from the normalization integral of the shape function over the momentum of the QQ in the quarkonium rest frame. This significantly enhances the predictive power of the shape function formalism, compared to existing studies where models for the shape functions had to be chosen arbitrarily. Second, the UV divergences induce momentum-dependent mixings between channels, which modify the matching conditions that are needed to determine the short-distance coefficients. This makes the matching conditions originally obtained in ref. [21] invalid beyond tree level, because there the mixing effects were not taken into account. Since current-day phenomenological studies of quarkonium production are usually carried out at next-toleading order (NLO) in α s , a correct derivation of the matching conditions incorporating mixing effects is necessary.
In this work, we consider the inclusive production of P -wave heavy quarkonia in the shape function formalism at NLO accuracy. We focus on the simpler case of production of P -wave quarkonia, which involves only one unknown nonperturbative matrix element, and the renormalization of the matrix elements is less complicated compared to the S-wave case. The P -wave production rate is also important in understanding the feeddown contributions in S-wave quarkonium production. In order to correctly incorporate the renormalization of NRQCD matrix elements, we compute the shape functions in perturbative QCD at NLO accuracy, which lets us obtain one-loop matching conditions for the shape function formalism that are necessary for computing the short-distance coefficients. This also gives us the large-momentum asymptotic forms of the shape functions, which, together with the normalization condition, help strongly constrain the nonperturbative form of the shape function and significantly reduce the model dependence. Furthermore, we can compute the shape functions in terms of quarkonium wavefunctions and universal gluonic operator vacuum expectation values using the potential NRQCD (pNRQCD) formalism developed in refs. [14,16,[23][24][25][26][27][28][29]. Although such a calculation does not yet lead to a first-principles determination of the shape functions due to our lack of knowledge of the nonperturbative dynamics of gluons, the universality of the gluon operator vacuum expectation values can enhance the predictive power of the formalism, similarly to the pNRQCD calculations of NRQCD matrix elements. Based on the results for the shape functions and the corresponding short-distance coefficients we obtain in this work, we resum the kinematical corrections to inclusive χ c and χ b hadroproduction rates at small and large transverse momentum.
This paper is organized as follows. In sec. 2 we introduce the NRQCD and shape function formalisms for P -wave quarkonium production. In sec. 3 we compute the shape functions in perturbative QCD, which are necessary in deriving loop-level matching conditions for the shape function formalism. In sec. 4 we compute the shape functions nonperturbatively in pNRQCD and obtain expressions in terms of quarkonium wavefunctions and universal gluonic operator vacuum expectation values. In sec. 5 we establish the matching conditions for the shape function formalism that allow us to obtain short-distance coefficients from the known results in NRQCD. By using these results we discuss phenomenological applications of the shape function formalism for hadroproduction of χ c and χ b at large and small transverse momentum in sec. 6. We conclude in sec. 7.

Shape function formalism for P -wave production
We first review the NRQCD factorization formalism for production of a χ QJ for J = 0, 1, and 2, where Q = c or b. At leading order in v, the inclusive production cross section of χ QJ is given by [1] σ[χ QJ (P )] = (2J + 1) c 3 P [1] where P is the momentum of the χ QJ , c 3 P 0 ) and O χ Q0 ( 3 S [8] 1 ) are NRQCD long-distance matrix elements that describe the evolution of a QQ in a specific color and angular momentum state into χ Q0 +anything. We have used the heavy-quark spin symmetry relations, which are valid up to corrections of order v 2 , to write the χ QJ matrix elements in terms of the ones for χ Q0 . The definitions of the NRQCD matrix elements read where d is the number of spacetime dimensions, ψ and χ are Pauli spinor fields that destroy and create a heavy quark and antiquark, respectively, D = ∇ − igA is the gauge-covariant derivative, A is the gluon field, ← → D is defined through the relation ψ † ← → D χ = ψ † Dχ−(Dψ) † χ, σ is a Pauli matrix, P χ Q0 is a projection onto states that include the quarkonium χ Q0 at rest, and T a is a SU (3) color matrix. In the color-octet operator O χ Q0 ( 3 S [8] 1 ), the adjoint Wilson line Φ (x) = P exp[−ig ∞ 0 dz · A(x + z)], where P is the path ordering for color matrices and is an arbitrary lightlike direction, is inserted to ensure the gauge invariance of the matrix element [30][31][32]. The QQ bilinears on the right and left of the projection operator P χ Q0 create and destroy a QQ in a specific color and angular momentum state. Equation (2.1) follows from a factorization conjecture, where the nonperturbative longdistance physics of scales below the heavy-quark mass m is encoded in the matrix elements, while the short-distance coefficients depend only on the short-distance process of the perturbative production of QQ. In this case, a QQ in either a color-singlet or color-octet state can evolve into a color-singlet quarkonium by emitting an arbitrary number of light particles with soft momentum of scales mv and lower. Hence, the QQ produced in a shortdistance process can have a relative momentum of order mv with respect to the produced quarkonium. The nonvanishing of this relative momentum can have a significant effect on quarkonium production observables, if the QQ is produced in the short-distance process predominantly near the boundary of phase space. In this case, logarithmic corrections to the short-distance coefficients from radiation of soft gluons can become enhanced, which must be matched with the effect of soft gluon emission in the NRQCD matrix elements. Also, in experimental observables kinematical cuts on the energy or momentum of the quarkonium are often taken, which can be sensitive to the order-mv change in the quarkonium momentum. In the NRQCD factorization formula in its usual form, the effect of the relative momentum between the QQ and the quarkonium is not taken into account, because the scales mv and lower do not affect the short-distance coefficients. Hence, in NRQCD, the momentum of the QQ produced in the hard process is identified with the momentum of the quarkonium. Therefore, predictions based on the factorization formula in eq. (2.1) can suffer from large corrections from kinematical effects associated with the motion of the QQ in the quarkonium rest frame.
The kinematical effects can be taken into account by including contributions from certain higher dimensional operators that involve total derivatives of the QQ bilinears in the factorization formula. The motion of the QQ along a momentum l gives rise to the following forms of higher dimensional matrix elements [18][19][20]: where the derivatives D act to the right and n is a positive integer. The operators involving derivatives on the left of the projector P χ Q0 can be obtained by using Hermitian conjugation and the invariance of the vacuum. Note that matrix elements of this kind do not appear in exclusive production due to conservation of energy and momentum. As have been shown in refs. [21,22], near the boundary of phase space, contributions from the above matrix elements associated with a lightlike momentum l collinear to the quarkonium momentum P become enhanced. These matrix elements can be obtained from moments of the "shape functions" defined by where δ(x) is the Dirac delta function, and the + direction is defined along the quarkonium momentum P . We take the convention l ± = l 0 ± l z , where l z = P · l/|P | in the frame where the quarkonium three-momentum P is nonzero. From the above definitions, it is evident that they are formally normalized by 1 ) . The shape functions defined in eqs. (2.4) can be interpreted as the probabilities for a QQ to evolve into a quarkonium after emitting soft particles with total momentum l. The requirement that NRQCD factorization holds after inclusion of the higher dimensional matrix elements in eqs. (2.3) constrains the shape functions to be defined only for l + > 0. This is because a negative l + , corresponding to the case where the QQ absorbs energy before evolving into a quarkonium, implies that the soft interactions between the QQ and the environment are not disentangled, and factorization is explicitly broken [33].
While the lightlike direction of the gauge-completion Wilson line is arbitrarily chosen in the definition of the color-octet matrix element, its origin is the direction of the heavy quarkonium momentum in a boosted frame such as the hadron CM frame [30,31], just like the direction of the lightlike momentum l. Although the gauge-completion Wilson line will not play a role in the phenomenological analysis of χ QJ production at NLO accuracy, because diagrams that involve the gauge-completion Wilson line begin to appear from two loops [30][31][32]34], for definiteness we will take the direction to be the same as l.
The factorization assumption leads to the following form of the factorization formula for the shape function formalism where the s N (P + l) are the short-distance coefficients in the shape function formalism, which must be functions of the QQ momentum P + l. This follows from the fact that the shape function formalism is obtained by resumming contributions from higher dimensional NRQCD matrix elements of the form given in eqs. (2.3), which implies that the coefficients s N (P + l) can be obtained from the standard NRQCD matching procedure by using the perturbative cross sections σ[QQ(P +l)], formally by expanding in powers of l. In ref. [21], it has thus been suggested that the short-distance coefficients in the shape function formalism are simply given by s N (P + l) = c N (P + l). However, this can become invalid beyond tree level, because UV divergences in the NRQCD matrix elements induce mixing between the two channels, which occur in different forms in NRQCD and shape function formalisms.
Let us compare the matching conditions for the two formalisms. The matching conditions that determine the NRQCD short-distance coefficients read (P ), which cancels the scale dependence of the coloroctet matrix element in the factorization formula. This kind of mixing occurs also in the shape function formalism, but in a different, l-dependent manner. The matching conditions that determine the short-distance coefficients s N (P + l) in the shape function formalism read (l + ) at loop level. We will show later that to one-loop accuracy, we have S 1 ) , because a QQ can evolve into a QQ in the same color and angular momentum state without emitting any gluons, but the same process cannot occur by exchange of a single gluon due to conservation of color. These relations give (P ) for the color-octet channel. The s 3 P [1] J (P ) for the color-singlet channel on the other hand, involves the color-octet shape function S (l + ) acquires at one loop a contribution that is proportional 0 ) . As we will see in sec. 3, an explicit calculation gives (P + l) with eq. (2.8). In order to properly regularize the divergences associated with the integral over l + in the factorization formula, the d-dimensional expression must be used in the matching condition, which requires an explicit calculation. Explicit calculations of the color-octet shape function S 1 ) , will be presented in section 3, so that we can obtain the short-distance coefficients s N from the NRQCD short-distance coefficients c N by comparing eqs. (2.6) and (2.7).
The result for the color-octet shape function in eq. (2.8) in d = 4 dimensions can be understood without explicit calculations in the following way. The normalization condition must reproduce the UV divergence in the color-octet matrix element when integrated over l + . The UV divergence in the color-octet matrix element gives the following evolution equation [1,29] where Λ is the scale for the renormalized matrix element O χ Q0 ( 3 S [8] 1 ) . In the perturbative calculation of this anomalous dimension, the divergence occurs from the integral over the momentum of the gluon emitted by the color-octet QQ as it evolves into a color-singlet Pwave state. Hence, the UV divergence in the normalization integral  (l + ) must come from the l + dependence of the shape function 1 . This implies that at large l + mv 1 In the case of the color-singlet matrix element, its two-loop UV divergence comes from renormalization of the derivative of the P -wave wavefunction at the origin [35,36], and so is unrelated to the l+ dependence of the color-singlet shape function.
the nonperturbative color-octet shape function must take the following asymptotic form so that the integral depends only on one scale l + , the above expression for the color-octet shape function is valid at d = 4 for all l + when the χ Q0 state is replaced by the perturbative QQ( 3 P [1] 0 ) state; hence we obtain eq. (2.8).
With the short-distance coefficients obtained from perturbative calculations, P -wave quarkonium cross sections can be computed in the shape function formalism once the nonperturbative shape functions S (l + ) are determined. In the case of the color-singlet shape function, the fact that effects of emissions of order mv gluons are suppressed by v 2 at the amplitude level leads to the observation that S 0 ) up to corrections suppressed by v 4 (the same conclusion can be obtained from the vacuum-saturation approximation, as was done in ref. [21]). On the other hand, it has been argued that the nonperturbative color-octet shape function S cannot be computed except within models. Even so, the model dependence can be significantly reduced by using the large-l + asymptotic form (2.10) and the normalization condition 1 ) . As the normalization condition requires the integral to be infrared (IR) finite, the integrand must not grow like 1/l + as l + → 0, while it must reproduce the asymptotic form for large l + . Although describing the behavior of the coloroctet shape function at small l + may still require models, model parameters can be fixed from the color-octet NRQCD matrix element by using the normalization condition. As renormalization of NRQCD matrix elements is usually done in the MS scheme, we have where we define the MS scheme at scale Λ by subtracting the UV pole and rescaling the MS scale µ through the relation µ 2 = Λ 2 e γ E /(4π). This expression is not very useful as it is, because it requires use of a model shape function in d dimensions. Instead, we may regulate the UV divergence by cutting off the l + integral as which is valid to order 0 . For large enough l max + , the last term can be computed in perturbation theory, because the large-l + behavior is fixed by the renormalization of the color-octet matrix element. That is, (l + ).

If we define
(2.14) where the last term subtracts the UV divergence in the l + integral, the normalization condition for the shape function can be written as where now every term is finite at d = 4 and the dependence can be dropped. The term ∆(Λ, l max + ), which translates cutoff regularization to the MS scheme, can be computed perturbatively as a series in α s . At leading nonvanishing order, it satisfies (l + ) at both large and small l + . The predictive power of the shape function formalism can be further improved by computing the nonperturbative shape functions in pNRQCD, similarly to what have been done for the NRQCD matrix elements, to express them in terms of products of quarkonium wavefunctions and vacuum expectation values of gluonic operators [14,16,28,29]. We will compute the shape functions in the pNRQCD formalism in section 4.
Before concluding this section we need to discuss the transformation of the shape function under Lorentz boosts. This is necessary because while the NRQCD operator definitions of the shape functions require them to be computed in the rest frame of the heavy quarkonium, similarly to the NRQCD matrix elements, the short-distance coefficients often need to be computed in a frame where the quarkonium has energetic spatial momentum. It is clear that from the definition of the shape functions, dl + S χ Q0 N (l + ) is Lorentz invariant under boosts in the + direction. Hence, the transformation of the shape function itself can be inferred from the boost property of l + . Under boosts from the quarkonium rest frame, l + transforms like where l * is the momentum l in the quarkonium rest frame, and l + on the left-hand side is the + component of l in the boosted frame; likewise, P * is the quarkonium momentum in the quarkonium rest frame, and P + is the + component of P in the boosted frame. Note that in the quarkonium rest frame, P * + = √ P 2 , which coincides with the quarkonium mass.
Throughout this paper, we will refer to the momentum l + in the quarkonium rest frame as l * + . Due to the invariance of dl + S χ Q0 N (l + ) under boosts, the first term of the left-hand side of eq. (2.15) is also invariant, as long as l max + boosts like l + . We will also see in the calculation of the scheme conversion ∆(Λ, l max + ) that the dependence on the boost cancels between the explicit dependence on l max + and the boost dependence coming from the scaling violation of the perturbative shape function S

Shape functions in perturbative NRQCD
In this section we compute the shape functions S (l + ) in perturbation theory, which are necessary for obtaining the matching conditions that determine the short-distance coefficients s N of the shape function formalism. As we have stated in the previous section, these objects are defined in the rest frame of the final-state QQ, and will be computed in this section in that frame as functions of l * + . The calculation of these shape functions are done in the same way as NRQCD matrix elements, except that the + component of the relative momentum of the QQ produced and annihilated on the left and right of the projection operator is constrained to be l * + . At order α 0 s , these shape functions are given by the squared amplitudes for QQ → QQ, where the color and angular momentum states of the initial and final state QQ are constrained to be 3 P 1 . Since they vanish unless l * = 0 and the color and angular momenta of the initial and final states are same, we obtain Now we consider the order-α s contributions. At order α s , the shape functions are given by the squared amplitudes for QQ → QQ with a gluon exchange between any of the heavy quark and antiquark lines. The Feynman diagrams for the order-α s contributions are shown in fig. 1. Note that the contribution from the gluon field in ← → D cancels in the operator definition of the color-singlet shape function, so that there are no diagrams involving direct gluon attachments to the QQ bilinear operator. While same diagrams appear in the order-α s calculation of shape functions and NRQCD matrix elements, contributions from the diagrams where the gluon crosses the final-state cut, like the first diagram in fig. 1, differ between the NRQCD matrix elements and the shape functions, because in the shape functions the momentum component l * + is not integrated over. If both the initial and the final-state QQ are in color-singlet states, such diagrams vanish due to conservation of color. In the case of S QQ( 3 S   vanish because they involve symmetric linear combinations of the totally antisymmetric SU (3) structure constants. The contributions from the remaining diagrams are same for both the shape function and the NRQCD matrix element. Hence, the tree-level results for In the case of S , the contribution from the gluon attachment on the heavy quark line cancels the contribution from the gluon attachment on the heavy antiquark line on the same side of the cut. The remaining diagrams cannot produce a color-octet QQ in the final state, and so they do not contribute to S holds through order α s , and for the same reason, 0 ) vanishes. While they can become nonzero from corrections of higher orders in α s , they will be suppressed by powers of v, because they involve dynamical gluons of order mv crossing the cut. Hence, they vanish at current order in v.
The only nontrivial contribution comes from S , we only need to consider the diagrams where a gluon crosses the cut, because otherwise the color-octet operator cannot produce a color-singlet QQ in the final state. We set the momenta of the quark and antiquark in the final state to be (E, q) and (E, −q), respectively, with E = q 2 /(2m). We use the NRQCD Feynman rules in ref. [37]. The contributions from gluon attachments to the Q and theQ on one side of the cut as shown in fig. 2 involve the factor ig(k + 2q)T a 2m where k is the momentum of the gluon, and we neglect the terms q 2 /(2m), (q + k) 2 /(2m) and (q − k) 2 /(2m) compared to k 0 in the heavy quark and antiquark propagators, because they are suppressed by a power of v. The product of the factor q from the gluon vertices and the Pauli matrix σ coming from the operator definition of the shape function can be decomposed into the trace, antisymmetric, and traceless symmetric parts as which is valid in d dimensions. In 3 spatial dimensions the three terms correspond to the irreducible angular momentum tensors for total angular momentum 0, 1, and 2. If the QQ in the final state is in the 3 P [1] 0 state, only the first term survives, and the factor q · σ can be identified as the tree-level matrix element of the QQ bilinear operator between the vacuum and the QQ( 3 P [1] 0 ) state. Then, the shape function is given by where Λ is the MS scale. Here, the 2πδ(k 2 ) comes from the gluon crossing the cut, and the δ(l * + − k + ) comes from the delta function in the definition of the shape function. The factor d − 2 comes from the sum over the final-state gluon polarizations, and the factor 1/(d − 1) comes from the projection onto the 3 P where C is given by Expanding C in powers of to linear order yields C = 1 + log is the result for the color-octet shape function S (l * + ) valid to order α s . Note that we recover the known result for the unrenormalized (bare) color-octet matrix 1 ) bare at order α s by integrating over l * + : This result also lets us compute the scheme conversion ∆(Λ, l max + ). From the definition given in eq. (2.14), we have at order-α s accuracy which is our result for the scheme conversion in the quarkonium rest frame. Note that this result reproduces the evolution equation in eq. (2.16). The scheme conversion vanishes at this order if we set l * max 1 ) at scale Λ = l * max + e 1/6 . From e −1/6 ≈ 0.85, we see that the cutoff l * max + is numerically close to the MS scale. In phenomenological studies of heavy quarkonium production, the MS scale Λ for the color-octet matrix element is usually chosen to be the heavy quark mass m, and in this case the corresponding cutoff l * max + equals m × e −1/6 . Consistently with the assumption that perturbative matching calculations are valid at a suitably chosen factorization scale Λ, we will also assume that the behavior of the shape function is perturbative for l * + at around or above l * max + , so that the l * + dependence of the scheme conversion and the shape function can be described by the perturbative calculation in this section.
We can now explicitly address the boost dependence of ∆(Λ, l max + ). In a boosted frame, eq. (3.10) becomes which we obtain by using l + = (P + /P * + )l * + . The factor (P + /P * + ) 2 in the integrand is due to the scaling violation of S , the result is exactly same as ∆(Λ, l * max + ) computed in the rest frame. Now that we have the perturbative results for the shape functions valid to orderα s accuracy, the short-distance coefficients s N for the shape function formalism can be obtained at one-loop level from the matching conditions in eqs. (2.7) by plugging in the results for the shape functions in eqs. (3.2), (3.3), and (3.7). Expressions for the s N in terms of the NRQCD short-distance coefficients c N will be given in section 5.

Nonperturbative analysis of shape functions in pNRQCD
The pNRQCD formalism [23,24,27] has been proven useful in analyses of inclusive quarkonium production [14,16,28,29]. This formalism provides expressions for NRQCD matrix elements as products of quarkonium wavefunctions and vacuum expectation values of gluonic operators, and reproduces the known results in terms of quarkonium wavefunctions in the case of color-singlet matrix elements. The vacuum expectation values of gluonic operators, often called gluonic correlators, are defined in terms of Wilson lines and gluon field-strength tensors. Although first-principles calculations of gluonic correlators appearing in NRQCD matrix elements have not yet been done, the fact that they are universal, that is, they are independent of the specific quarkonium state including the heavy quark flavor and radial excitation, significantly enhances the predictive power of the nonrelativistic effective field theory formalism. We refer the readers to refs. [14,16,28,29] for detailed discussions of pNRQCD calculations of NRQCD matrix elements and phenomenological applications in inclusive quarkonium production.
We can expect to obtain similar results from pNRQCD calculations of shape functions: shape functions can be written in terms of quarkonium wavefunctions and l + -dependent gluonic correlators. If the l + dependences are solely contained in the gluonic correlators, then we can expect the l + dependences of the shape functions to be universal among shape functions of quarkonium states that differ by radial excitation or heavy quark flavor. As the model dependence of the shape function formalism comes from the l + dependence of the nonperturbative shape functions, this would greatly enhance the predictive power of the shape function formalism.
In this section, we compute the shape functions S (l + ) in the pN-RQCD formalism for P -wave quarkonium production developed in refs. [28,29]. Similarly to the calculations of NRQCD matrix elements, we will compute them in the quarkonium rest frame as functions of l * + . We derive the results nonperturbatively following the development of the formalism in refs. [28,29], which results in expansions in powers of v and Λ QCD /m, while the dynamics of the gluon and light quarks are kept nonperturbative. As have been done in refs. [28,29] we assume that the χ Q states are strongly coupled, so that the scale mv 2 is smaller than the typical energy gap of gluonic excitations. We will discuss the possibility of extending to the weak coupling case after obtaining the results in the strongly coupled case.
We first briefly review the calculation of NRQCD matrix elements in the formalism developed in refs. [28,29]. In the pNRQCD formalism, the matrix elements are given by is the wavefunction of the quarkonium Q at leading order in v, x 1 and x 2 are the positions of the heavy quark and antiquark, respectively, r = x 1 − x 2 , and R = 1 2 (x 1 +x 2 ). Similar relations hold for the primed coordinates.
is a contact term that is determined by matching NRQCD and pNRQCD. The contact terms for the NRQCD matrix elements in eq. (2.1) are given by where we discard contributions that vanish when applied to wavefunctions in P -wave states.
The tensor E ij 11 is defined by where T andT indicate time and anti-time orderings of the operators, E a,i (t, x) is the chromoelectric field at time t and spatial position x, and Φ 0 (t 1 , x; t 2 , x) is an adjoint Wilson line in the temporal direction defined by The configuration of the Wilson lines in the definition of the tensor E ij 11 is shown in fig. 3. Similarly to eq. (3.5), we can decompose the product σ k ∇ i r into the trace, antisymmetric, and traceless symmetric parts; if we apply the contact term to a wavefunction in the 3 P 0 state, the antisymmetric and traceless symmetric parts vanish. Hence, we have where E 11 is the trace of the tensor E ij 11 , which can be written as Here we suppress the spatial positions of the operators as they are all evaluated at x = 0.
The relation between E 11 and the quantity E defined in refs. [28,29] is given by E 11 =  E in d spacetime dimensions 2 . By plugging in the results for the contact terms to the pNRQCD formula in eq. (4.1), we obtain where R(r) is the radial wavefunction for the χ Q0 state. While the nonperturbative value for E 11 has not yet been computed from first principles, the fact that it is a universal quantity that does not depend on the radial excitation or the heavy quark flavor enhances the predictive power of the nonrelativistic effective field theory formalism. For instance, we have for any 3 P 0 quarkonia. Note that E 11 has a logarithmic scale dependence that reproduces the evolution of the color-octet matrix element [28,29].
We note that the derivatives ∇ r on the left and right of the δ (d−1) (r) in the contact terms in eqs. (4.2) come from the QQ bilinears on the left and right of the projection operator. Similarly, the two chromoelectric fields in E 11 also come from the QQ bilinears on the left and right of the projection operator. By using these points, it is straightforward to obtain pNRQCD expressions for the shape functions. Rather than working directly with the definitions in eq. (2.4), we define the Fourier transforms where now the QQ bilinears on the right of the projection operator are displaced by b, where b lies in the − direction. The shape functions in momentum space can be recovered from S Other than that the QQ bilinears on the left and right of the projection operator are computed at different spacetime positions, the pNRQCD calculations of the shape functions are done in the same way as NRQCD matrix elements. We havẽ where the contact terms for the shape functions acquire dependence on b − . The contact terms that survive when acting on wavefunctions in the 3 P 0 state read Compared to the expressions for the contact terms for NRQCD matrix elements, the only difference in the shape function formalism is that inẼ / 11 (b − ), the chromoelectric field on the right and the Wilson lines associated with it are displaced by b. The displacement does not affect the derivatives ∇ r , because they only depend on the relative coordinates between the Q andQ, and they apply to eigenstates of NRQCD. Since the contact term for the color-singlet shape function is independent of b − and is the same as the contact term for the color-singlet NRQCD matrix element, we havẽ On the other hand, the Fourier transformed color-octet shape function acquires dependence (4.15) Then the shape function can be written as where Similarly to the case of NRQCD matrix elements, the quantity E / 11 (l + ) is universal, and do not depend on the radial excitation or the heavy quark flavor of the quarkonium state. In particular, we have which is valid for any 3 P 0 quarkonia. This, together with heavy quark spin symmetry, implies that the l + dependence of the color octet shape function is universal for all P -wave heavy quarkonia, including all radial excitations of χ c and χ b . This significantly reduces the model dependence in the phenomenological analysis of χ c and χ b production in the shape function formalism. An important consistency check of the pNRQCD calculation is to compare the result with the perturbative calculation. If we compute E / 11 (l * + ) in perturbative QCD at leading nonvanishing order in α s , we obtain where C is defined in eq. (3.8). By plugging in this expression into the pNRQCD expression for the shape function (4.16) and using O χ Q0 ( 3 P which exactly reproduces the order-α s calculation of S demonstrating the validity of the pNRQCD result at one-loop level. This is also consistent with the result for the scheme conversion ∆(Λ, l * max As have been discussed in section 3, the UV divergence on the left-hand side is regulated by the cutoff l * max + , so that the nonperturbative shape function S (l * + ) can be computed in d = 4 dimensions. This allows us to use the UV-regulated normalization condition to constrain the model dependence of the nonperturbative shape function.
So far the pNRQCD calculations of the NRQCD matrix elements and shape functions presented in this section were based on the assumption that the χ Q0 state is strongly coupled, so that the ultrasoft scale, where the gluons have energies and momenta of order mv 2 , does not play a role. In the weakly coupled case, the chromoelectric dipole interaction will mix color-singlet and color-octet QQ states through exchange of ultrasoft gluons, so that the contribution from the color-octet wavefunction must be included. As have been discussed in refs. [38,39], this interaction is suppressed by v 3/2 compared to the kinetic terms in the pNRQCD Lagrangian, so that for heavy quarkonium states, the contribution from the color-octet wavefunction is suppressed by v 3 compared to the color-singlet part. The color-octet wavefunction can appear in the color-octet NRQCD matrix element and the shape function through the contact terms at leading order in 1/m, which is enhanced by 1/v 2 compared to the contact terms that act on the color-singlet wavefunction [eqs. (4.2) and (4.11)]. As a result, in the weak coupling case, the extra contribution coming from the color-octet wavefunction is suppressed by at least a power of v compared to the pNRQCD expressions for NRQCD matrix elements and shape functions in the strongly coupled case. Hence, we expect the pNRQCD results obtained in this section to also hold in the weakly coupled case at leading nonvanishing order in v.

One-loop matching conditions
Based on the calculations of the perturbative shape functions, we now derive the matching conditions that allow determination of the short-distance coefficients s N for the shape function formalism.
As have been discussed in section 2, the matching coefficients are computed from the cross sections for inclusive production of QQ in specific color and angular momentum states. The cross section for production of QQ in the 3 S where the first line is the NRQCD expression divided by 2J +1, and the last line follows from the result for the perturbative shape function in eq. In principle, we can obtain expressions for s 3 P 1 ) contains an IR pole, and hence, the NRQCD side of the expression depends on the ordercontribution to the short-distance coefficient c 3 S [8] 1 (P ) that is usually not available. The same IR pole that matches the one on the NRQCD side of the expression does appear on the shape function formalism side, which arises from the integral in the last line of eq. (5.2) in the region l + → 0. Hence, the terms involving IR poles on both sides can be matched exactly if we rewrite this integral as where the unrenormalized ("bare") matrix element is given by eq. (3.9). Note that the integral over l + on the right-hand side is IR finite because the integrand factor in the parenthesis vanishes as l + → 0, and the only IR pole is isolated in the perturbative color-octet matrix element. This expression is still not completely satisfactory, because NRQCD matching calculations are usually done in the MS scheme, while the above expression involves an unrenormalized matrix element. Rather than subtracting the pole to carry out the renormalization in the MS scheme, we can use the result for ∆(Λ, l max + ) we obtained earlier in eq. (3.10) to cut off the UV-divergent integral. By using eqs. (2.15) and (3.10), we have at one-loop level where θ(x) is the step function that vanishes for x < 0 and equals 1 for x ≥ 0, and Λ is the 1 ) . By using this we can rewrite the color-octet shape-function term in the matching condition as where now every term is UV finite, and the IR pole is isolated in the MS-renormalized matrix element at the scale Λ = l * max + e 1/6 . By plugging this expression to eq. (5.2) we obtain If we solve this for s 3 P [1] J (P ), the terms involving IR poles from the color-octet matrix element cancel exactly, so everything can be computed at d = 4. By using the result for (l + ) in eq. (3.7), we obtain where c singular (5.8b) The labels "singular" and "pole" will be explained later. Since the poles cancel in the last line of eq. 1 ), which then evolves into a P -wave color-singlet QQ by emitting a soft gluon. The soft gluon emission amplitude from a heavy quark line with momentum p 1 is given bȳ where we retained only the leading contribution in the limit where k is soft, and in the last equality we anticommuted p / 1 to the left and used the equation of motion. Similarly, the infrared divergent contribution coming from the heavy antiquark line with momentum p 2 is given by Hence, the soft-gluon emission gives rise to an infrared divergent contribution that is given by multiplied on both sides of the cut, and integrated over the phase space of the soft gluon.
We can see that this soft gluon factor is proportional to 1/|k|, so that the phase space integral involves the singular integral d d−1 k |k| −3 . After a straightforward calculation, at leading nonvanishing order in expansion in powers of q = (p 1 − p 2 )/2 compared to p 1 + p 2 , the singular contribution to the QQ( 3 P

[1]
J ) cross section is given by which, aside from an order-0 finite piece, matches exactly the contribution from the c singular 3 P The soft-gluon emission that is subtracted from the color-singlet short-distance coefficient instead appears in the color-octet shape function contribution, where the nonperturbative effect of the soft gluon emission is encoded in the shape function. To see this explicitly, we can rewrite the color-octet shape function term in the cross section formula (5.9) as ( 5.15) We see that the second term in eq. (5.7) is equivalent to the last line of the above expression, with an opposite sign and the shape function replaced by its large-l + asymptotic behavior given in eq. (2.10). By using this we can rewrite eq. (5.9) as Here, the first line is just the cross section in NRQCD, while the remaining terms arise from the deviation of the nonperturbative shape function from its asymptotic form. That is, the integral over l + in eq. (5.16) encode the nonperturbative corrections to the NRQCD factorization formalism that arise from resummation of kinematical effects from the motion of QQ relative to the quarkonium. Although eqs. (5.9) and (5.16) are in principle equivalent expressions for the χ QJ production cross section in the shape function formalism, they have different phenomenological implications, as we will now briefly describe. In order to compute cross sections from eq. (5.9), the short-distance coefficient s 3 P at tree level. However, since usually in NRQCD calculations both the color-singlet and color-octet short-distance coefficients are computed to same accuracy in α s , there is always a mismatch of the orders in α s of soft-gluon effects between the color-singlet and color-octet channels, due to the one-loop correction included in c 3 S [8] 1 . This mismatch can have a significant effect, as the soft gluon emission involves IR divergences which make the cross sections singular near threshold. If we use eq. (5.9) to compute χ QJ cross sections, this singular contribution is subtracted from the 3 P

[1]
J shortdistance coefficient, and added back to the second term in the square brackets of eq. (5.9) in the form of eq. (5.15). In this way, the singular soft-gluon emission effects are contained entirely in the color-octet shape function contribution in the form given in eq. (5.15), where the short-distance coefficient c 3 S [8] 1 can now be computed to same order in α s for every term. This allows the singular soft-gluon emission effects to be computed to consistent accuracy in α s . Note that, since the integral over l + in eq. (5.15) is infrared finite even when the nonperturbative shape function is replaced by its asymptotic form, this matching of singular soft-gluon effects can be done without knowledge of the nonperturbative form of S χ Q0

S
[8] 1 (l + ). We will investigate this point in more detail in section 6.1.
The form of the cross section given in eq. (5.16) is suitable for computing the nonperturbative corrections to the NRQCD factorization formalism from the nonperturbative shape function. We note that if we were to neglect the deviation of the nonperturbative shape function from its asymptotic form, the last line of eq. (5.16) vanishes and we recover the cross section in NRQCD factorization. This is consistent with the fact that, if the color-octet shape function is given by the asymptotic form in eq. (2.10), all of the higherdimensional matrix elements in eq. (2.3b) are proportional to scaleless power divergences ∞ 0 dl + l n−1 + for n ≥ 1, which vanish in dimensional regularization (in the case of higherdimensional color-singlet matrix elements, they vanish in the form ∞ 0 dl + l n + δ(l + ) = 0). The nonperturbative corrections to the p T -differential cross sections coming from the last two lines in eq. (5.16) will be computed in section 6.2 using models for the nonperturbative shape functions, which are constrained by the asymptotic form (2.10) and the normalization condition given by eq. (4.21).

Phenomenological applications
Based on the one-loop level matching conditions and the shape functions we have established in the previous sections, we compute χ c and χ b cross sections in the shape function formalism and compare the results with NRQCD predictions. As we have explained in sec. 1, understanding the χ c and χ b cross sections are not only important for comparing predictions with cross section measurements in experiment, but they are also important for correctly treating the feeddown contributions in S-wave production rates. We will show in this section that the shape function formalism can modify the transverse-momentum dependent cross sections of P -wave heavy quarkonia, which may have important implications in heavy quarkonium production phenomenology.

Matching of soft-gluon emission effects at large transverse momentum
We first investigate the phenomenological application of the cross section formula in eq. (5.9) for mixing of soft-gluon emission effects between color-singlet and color-octet channels. Since the effects of soft-gluon emission near threshold is most prominent at large p T , we consider the cross section in the fragmentation approximation, which describes the cross section at leading power in the expansion in powers of m 2 /p 2 T [40]. In this approximation, the cross section is given by where the sum is over gluon, light quarks and antiquarks,σ i (k) is the production rate of a parton i with momentum k, D i→χ QJ (z) is the fragmentation function, which is a function of the momentum fraction z = P + /k + , with the + direction defined along the quarkonium momentum. For notational convenience, we will write the parton cross sectionσ i as a function of z. The NRQCD factorization formula for the fragmentation function reads where the short-distance coefficients for the fragmentation functions can be computed perturbatively. The short-distance coefficients for the fragmentation functions begin at order α s , and the parton cross sectionsσ i begin at order α 2 s for proton-proton collisions. Hence, parton cross sections to order-α 3 s accuracy and fragmentation functions to order-α 2 s accuracy are needed to compute large-p T quarkonium cross sections at NLO in α s . The complete results for heavy quarkonium fragmentation functions to order-α 2 s accuracy can be found in refs. [41][42][43][44][45][46][47][48][49][50][51]. The quark fragmentation function d q→ 3 S [8] 1 begins at order α 2 s , while quark fragmentation into QQ( 3 P (z) to order-α 2 s accuracy read where µ f is the factorization scale for the fragmentation function, is the gluon splitting function [52][53][54][55]. The plus distribution is defined by The J-dependent constants Q J and functions P J (z) read (z) begins at order α s , and so the above expression contains corrections at NLO in α s . The term proportional to β 0 δ(1 − z) log µ f cancels the running of α s at lowest order, while the µ f dependence coming from the log µ f P gg (z) term cancels the µ f dependence in the gluon cross sectionσ g , so that the quarkonium cross section is independent of µ f . The distribution in the last line of eq. (6.3b) arises from soft-gluon emission, which is singular near the z = 1 threshold. Because the gluon cross sectionσ g (z) steeply rises with z, this singular term can have a significant effect on the cross section, especially at large p T . Note that such effect will not appear in d g→ 3 P (z) will then have more singular corrections at order α 3 s coming from multiple soft gluon emissions near threshold, and the corresponding correction in d g→ 3 P [1] J (z) will only appear at order α 4 s . While corrections from soft gluon emissions near threshold can be resummed [56], it is also important that they are treated consistently throughout the color-singlet and color-octet channels. This is especially the case in χ QJ production at large p T , where the contribution from the color-singlet channel can turn negative at values of p T much larger than the quarkonium mass, so that the cross section is given by the remnant of the cancellation between the color-singlet and color-octet channel contributions. Because of this cancellation, the cross section at large p T can become very sensitive to soft gluon emission effects near threshold, so that mismatch between the treatment of soft gluon emissions between color-singlet and color-octet channels can spoil the reliability of the perturbative expansion.
As we have argued in the previous section, the formulation of the cross section in the shape function formalism of the form given in eq. (5.9) can be used to treat the soft-gluon emission effects consistently between the color-singlet and color-octet channels. We will see explicitly by computing the fragmentation function in the shape function formalism that the singular parts of the fragmentation functions can be encoded entirely in terms of the color-octet short-distance coefficient. Let us write the fragmentation function in the shape function formalism in the form given in eq. (5.9) as D g→χ cJ +X (z) shape = (2J + 1) d s which we obtain from eq. (5.7). We used the boost invariance of l + /P + to write the integral in terms of l * + in the quarkonium rest frame. We need this expression at order-α 2 s accuracy, which comes from the order-α s contribution in d 3 S [8] 1 (z). Explicit calculation of the subtraction term is carried out as follows: where in the last line we used C F = 4/3 for N c = 3 and for any function f (z). From this we obtain We see that now the singular plus function contribution in d g→ 3 P [1] J (z) has disappeared, along with the dependence on the NRQCD factorization scale Λ. As we have shown in the previous section in eq. (5.15), this contribution reappears in the color-octet shape function term in eq. (6.6). If we rewrite eq. (6.6) using eq. (5.15) as we can now use the order-α 2 s expression for d g→ 3 S [8] 1 everywhere, so that the soft-gluon emission effects are consistently taken into account. Note that if we are only interested in the reorganization of the perturbative corrections arising from soft-gluon emission effects, we can replace the color-octet shape function in eq. (6.11) by its asymptotic form S 0 ) , because the integrand factor in the parenthesis vanishes as l + → 0, rendering the integral over l + IR finite. In this way, the contribution proportional to O χ Q0 ( 3 P [1] 0 ) given by is computed to order-α 2 s accuracy, can be interpreted as the color-singlet contribution to the fragmentation function augmented with order-α 3 s soft-gluon effects. Here, we have now set = 0, because the integral over l + is IR and UV finite.
To quantify the effect of the inclusion of the order-α 3 s soft-gluon effects to the colorsinglet contribution, we compute the color-singlet contributions to the cross sections in NRQCD and in the shape function formalisms defined by J (z) for NRQCD and D g→χ cJ +X (z)| singlet, shape is given in eq. (6.11) for the shape function formalism, and compare them with σ[χ QJ ] octet defined by (z) , (6.14) as functions of p T . Note that at this accuracy there is no quark fragmentation contribution to the color-singlet channel, because as we previously mentioned d q→ 3 P [1] J (z) begin at order α 3 s . In order to facilitate comparisons independently of the determinations of NRQCD matrix elements, we define the dimensionless ratios .

(6.15)
We compute the ratios R J for χ c1 and χ c2 production as functions of p T at central rapidity from proton-proton collisions at √ s = 13 TeV. We set the charm quark masses to be m = 1.5 GeV and choose the renormalization scale for the color-octet matrix element to be Λ = 1.5 GeV, as is usually done in phenomenological studies of charmonium production. We compute the gluon cross sectionσ g at NLO accuracy by using the code in ref. [57]. In order to make contact with existing NLO calculations of quarkonium cross sections in NRQCD, we neglect the order-α 5 s cross term that comes from the order-α 3 s gluon cross sectionσ g and the order-α 2 s color-octet fragmentation function. That is, we compute the g +O(α 4 s ) and d g→ 3 S [8] σ q can be computed to order-α 2 s accuracy because d q→ 3 S [8] 1 begins at order α 2 s . We computê σ g at NLO in α s with CTEQ6M parton distribution functions [58] at scale µ f = m T = p 2 T + 4m 2 , and compute α s at the same scale using the two-loop formula provided in ref. [58] with n f = 5 light quark flavors and Λ (5) QCD = 226 MeV. This calculation of the short-distance coefficients at leading power has been shown to agree well with the fixedorder calculation at large p T , and the fragmentation contributions completely dominate the cross section for p T above 100 GeV [11,59].
The results for R J=1 and R J=2 from NRQCD and the shape function formalism are shown in fig. 4. We concentrate on the J = 1 and J = 2 cases, because χ c0 is usually not measured in hadroproduction due to the tiny branching fraction for χ c0 → J/ψ + γ. Note that these ratios take negative values at large p T , because the color-singlet cross section involves a large negative contribution coming from the singular plus function term z/(1 − z) + in d g→ 3 P [1] J (z) associated with singlet-octet mixing induced by soft-gluon emission. To ensure the positivity of the χ cJ cross sections, the ratios must be greater than the ratio of matrix elements given by −m 2 O χ c0 ( 3 S 0 ) at Λ = 1.5 GeV determined from cross section measurements at hadron colliders is 0.043 at NLO accuracy [29,60]. We see that in the case of NRQCD, the ratios R J become more negative as p T increases and eventually make the cross sections turn negative, as they go below −0.043 when p T exceeds about 1.2 TeV for J = 1 and 0.9 TeV for J = 2. This problem does not occur until much larger p T in the shape function calculation, as the ratios R J are almost constant over a wide range of p T , and the positivity of the cross sections is ensured as the values of R J stay safely above −0.043 for much larger values of p T . As we have argued, this happens because the shape function formalism allows us to match the accuracy in α s of the soft-gluon emission effects between the colorsinglet and color-octet channels, so that the cancellation between the two channels can occur consistently without mismatch.
We note that although the shape function formalism eliminates the mismatch of softgluon emission effects between channels that mix under renormalization, this does not necessarily alleviate the need for resummation of threshold logarithms associated with soft-gluon emission, as current results are still based on fixed order-calculations. Even when carrying out the threshold resummation, the shape function formalism in the form of eq. (5.9) employed in this section will allow us to consider the resummation effects in the color-octet and color-singlet channels in a consistent manner, and ensure that there is no mismatch of soft-emission effects.

Nonperturbative corrections from shape functions
We now focus on the effect of nonperturbative corrections arising from the nonperturbative shape function that are neglected in usual applications of NRQCD factorization. That is, we make use of eq. (5.16) to compute the additional nonperturbative contribution coming from the deviation of the nonperturbative shape function from its asymptotic form. Because first-principles determination of the color-octet shape functions have not yet been done, this requires the use of model functions for the nonperturbative shape function. However, as we  have argued previously, the normalization condition (4.21) and the requirement that the nonperturbative shape function must coincide with the asymptotic form at large l + severely constrain the model dependence. Additionally, the pNRQCD results for the shape functions imply that the l + dependence of the shape functions must be same for any P -wave heavy quarkonia, independently of radial excitation or heavy quark flavor; this lets us use a single model function for all shape functions for production of any P -wave heavy quarkonium state.
We first attempt a crude estimate of the size of the nonperturbative corrections at large transverse momentum p T . We note that the factor S χ Q0 (l + ) asy in eq. (5.16) diminishes as l + increases, as the nonperturbative shape function must coincide with its asymptotic form when l + is in the perturbative regime. In the opposite case, when l + approaches zero, S χ Q0 3 S (l + ) asy diverges like 1/l + , while the nonperturbative shape function must either be finite or at least diverge slower than 1/l + in order to ensure the IR finiteness of its normalization. This implies that at very small values of l + , we can approximate S (l + ) asy . Hence, the additional nonperturbative contribution in eq. (5.16) can be estimated as (P ) as the extra momentum l reduces the available phase space. In the case of inclusive production at large p T , the p T -differential color-octet cross sections fall off like 1/p n T , with n ≈ 4. This, and the fact that at large p T the color-octet QQ is produced predominantly near threshold through gluon fragmentation leads to the following approximation for the l + dependence of the short-distance coefficient: where in the last equality we expanded in powers of l + to linear order compared to P + , because l NP + is much smaller than P + . By using this result, we can make the following estimate Note that the boost-invariant ratio l NP + /P + must be much smaller than 1/2, because P * + = 2m gives l max + /P + = e −1/6 /2 when Λ = m. If we set l NP + /P + = λ/2, with λ 1, the additional nonperturbative contribution to the cross section from the P -wave quarkonium cross section is estimated to be Compare this with the color-octet contribution to the cross section in the NRQCD factorization formalism given by where we used the pNRQCD result for the color-octet matrix element to write the cross section in terms of the color-singlet matrix element and the universal quantity E 11 . In the case of χ cJ , best-fit values of P -wave quarkonium matrix elements give E 11 = 1.17 at scale Λ = 1.5 GeV [29,60]. By setting α s = 0.25 at the scale of the charmonium mass, λ = 0.5, and n = 4, we find that σ[χ QJ (P )] NP is of roughly the same order of magnitude as σ[χ QJ (P )] octet in the case of charmonium production. In the case of bottomonium, σ[χ QJ (P )] NP is smaller than σ[χ QJ (P )] octet , because E 11 becomes larger at the scale of the bottom quark mass. Even though the additional nonperturbative effects from nonperturbative shape functions estimated in eq. (6.20) can be as large as order one compared to NRQCD calculations of the cross section, this by itself is inconsequential in large-p T quarkonium production phenomenology, because the additional contributions do not alter the p T dependence of the cross section and only brings in changes in the overall normalization. That is, the additional nonperturbative contributions can be absorbed into a redefinition of the color-octet matrix element. Because in phenomenological analyses the color-octet matrix elements are obtained from measured cross sections, and not computed from first principles, this redefinition does not affect the large-p T phenomenology in any way. It should be noted, however, that this estimate is valid only for the large-p T asymptotic behavior of the cross section.
At smaller values of p T , the shapes of the color-octet cross sections change, which would invalidate the assumptions that lead us to the estimate in eq. (6.20). Hence, at smaller values of p T , the additional nonperturbative contributions may change the p T shapes of the cross section significantly. In order to investigate these effects quantitatively, we need to compute σ[χ QJ (P )] NP explicitly with model shape functions.
Before we go on with the computation of the cross section, we discuss another nonperturbative effect associated with the heavy quarkonium mass that can be incorporated in the shape function formalism. As have been discussed in ref. [21], the shape function formalism allows matching the QQ phase space with the physical one involving the heavy quarkonium mass. In NRQCD factorization, the invariant mass of the QQ produced in the short-distance process can be different from the heavy quarkonium mass by order mv, due to the fact that the QQ can emit order-mv gluons before evolving into a quarkonium. In the shape function formalism, the kinematics of the emission of order-mv gluons is explicitly taken into account by the l + dependences in the shape functions and the corresponding short-distance coefficients. Hence, as have been suggested in ref. [21], in the calculation of the short-distance coefficients s N (P + l) and s N (P ), the QQ mass m QQ = √ P 2 can be chosen to be the heavy quarkonium mass, instead of m QQ = 2m. In principle, in the nonrelativistic power counting the difference between 2m and the heavy quarkonium mass is of order mv 2 . Because the shape functions describe the emission of order-mv gluons, and not necessarily order-mv 2 effects, the order-mv 2 effects can be consistently neglected from the short-distance coefficients by setting m QQ = 2m in the same way as done in NRQCD. However, because of the renormalon ambiguity in the heavy quark pole mass, the difference between the heavy quarkonium mass and 2m can exceed mv 2 and become numerically significant depending on the choice of the pole mass m. Setting m QQ to be the quarkonium mass avoids this issue and ensures that the QQ phase space matches the physical one within order mv 2 .
We now compute the χ c and χ b cross sections in NRQCD and shape function formalisms. Since we are interested in the nonperturbative effects from the shape function formalism at values of p T lower than the fragmentation regime, we compute the cross sections for the LHCb kinematics from 7 TeV pp collisions at the LHC, where data for P -wave quarkonium production rates relative to 1S quarkonium production rates are available for p T ≥ 2 GeV for χ c [61] and p T ≥ 6 GeV for χ b [62], which go below the quarkonium masses. The available p T ranges are, however, above the heavy quark pole masses which we set to be m = 1.5 GeV for charm and m = 4.75 GeV for bottom. The cross sections are computed in NRQCD and shape function formalisms as 1 )] dp T θ(l max We compute the NRQCD short-distance coefficients at NLO in α s numerically by using the FDCHQHP package [63] with CTEQ6M parton distribution functions [58]. Note that dσ[χ QJ ]/dp T | shape involve the same short-distance coefficients as NRQCD, but in this case we compute them with the QQ mass set to be the physical χ QJ mass, while in the NRQCD cross section the QQ mass is 2m, as is usually done in phenomenological studies. We compute the short-distance coefficients that enter the nonperturbative correction term dσ[χ QJ ]/dp T | NP at leading order in α s from the known parton cross sections dσ[ij → QQ( 3 S 1 )]/dp T at tree level as 24) where f i/p (x) is the parton distribution function for finding a parton i with momentum fraction x in a proton, y is the rapidity of the QQ, and y is the rapidity of the recoil X. The factor 2p T x 1 x 2 is the Jacobian arising from the change of integration variables from x 1 , x 2 , andt to p T , y, and y . The parton cross sections dσ[ij → QQ( 3 S [8] 1 ) + X]/dt are available as functions of the partonic Mandelstam variablesŝ,û, andt, which can be written in terms of s, p T , y, y , and P 2 = m 2 QQ aŝ where m T = m 2 QQ + p 2 T , x 1 = (m T e y + p T e y )/ √ s, and x 2 = (m T e −y + p T e −y )/ √ s.
Analytical expressions for the tree-level 2 → 2 parton cross sections for gg, gq, and qq initial states can be found in ref. [64]. The calculation of dσ[QQ( 3 S [8] 1 )]/dp T P →P +l is carried out in a similar way, except now the partonic Mandelstam variables are given bŷ / √ s, and l T is the transverse component of l. Note that there is some ambiguity in choosing the direction of the momentum l, because boosts along the beam direction will change the directions of P and l differently. We choose the directions of l and P to coincide in the frame where |P | = p T and |l| = l T . In this frame, l + = 2|l| = 2l T , and the boost to this frame from the rest frame of the χ Q is given by 1 )]/dp T P →P +l can be computed in a similar form as eq. (6.24), where now the l-dependent partonic Mandelstam variables are given in eqs. (6.26), and the Jacobian now takes a more complicated l-dependent form.
In order to compute dσ[χ QJ ]/dp T NP , the parton cross sections dσ/dt for the color-octet channel must be known as functions of the partonic Mandelstam variables. The analytical expressions for gg, gq, and qq initial states are available at tree level (order α 3 s ) in ref. [64]. At NLO (order α 4 s ), publicly available results are given only as p T -differential cross sections convolved with parton distribution functions and integrated over rapidity ranges, so that they cannot be used directly in eq. (6.23) to compute dσ[χ QJ ]/dp T NP . Fortunately, in the kinematical ranges that we are interested in this section, the effect of the NLO correction is small for the color-octet channel: as can be seen in ref. [8], the NLO K factor for the 3 S 1 )]/dp T when we compute the nonperturbative corrections dσ[χ QJ ]/dp T NP , while we compute the short-distance coefficients at NLO accuracy everywhere else.
Now we need a model for the nonperturbative shape function. As we have stated earlier, the nonperturbative shape function S χ Q0 3 S (l + ) asy for l + l max + , and normalized by 1 ) MS . For the normalization integral to be IR finite, the nonperturbative shape function must either be finite, or at least diverge slower than 1/l + as l + → 0. We consider the following model functions In model 1, the model shape function vanishes linearly as l + → 0, while in model 2 the model function becomes a finite constant. The model 3 shape function coincides the asymptotic form for l + > a + , while it vanishes for l + < a + . In model 4, the model function diverges as l + → 0 like 1/l 1−α + , more slowly than the asymptotic form, to ensure the IR finiteness of the normalization integral. We use the value α = 0.8 so that the model function slowly diverges at l + = 0, and recovers the asymptotic form at large l + . The model parameters a + , which transforms under boosts like l + , can be determined from the normalization condition. By using the value α s = 0.33 at scale 1.5 GeV, and the best-fit value for χ cJ matrix elements giving m 2 O χ c0 ( 3 S 0 ) = 0.043 at Λ = 1.5 GeV [29,60], we obtain the model parameters at the rest frame of the χ Q0 given by a * + = 0.74, 1.7, 0.64, and 2.6 in units of GeV for models 1, 2, 3, and 4, respectively. The shapes in l * + of theses model functions for fixed values of the model parameters are shown in fig. 5 compared to the asymptotic form.
In order to investigate the model dependence of the cross section in the shape function formalism, we compute dσ[χ QJ ]/dp T NP from the four model shape functions. In order to assess the impact of the nonperturbative corrections, we compute the dimensionless ratios of the color-octet channel contributions in the NRQCD and shape function formalisms where the numerator and denominator are defined by dσ[χ QJ ] dp T octet, shape Since we are interested in the behavior of dσ[χ QJ ]/dp T compared to the color-octet channel contribution in the usual NRQCD factorization formalism, we compute the short-distance coefficients dσ[QQ( 3 S 1 )/dp T at tree level (order α 3 s ) for the calculation of r Q octet (we will use the full NLO results in the calculation of cross sections for comparison with measurements). As we have stated earlier, we set the QQ mass to be 2m for the cross section in NRQCD, while we set it to be the quarkonium mass in the shape function formalism. The ratios depend on the ratio of NRQCD matrix elements m 2 O χ Q0 ( 3 S 0 ) . For charmonium we set it to be the best-fit value 0.043 at Λ = 1.5 GeV. For bottomonium, we compute the ratio for χ b0 at scale 4.75 GeV by using the universality of the gluonic quantity E 11 and its evolution equation [29]. That is, where the left-hand side equals E 11 /[(d − 1)N 2 c ] at scale 4.75 GeV, and the last line comes from the evolution of this quantity from scale 1.5 GeV to 4.75 GeV. Numerically, we obtain similar numbers if we use the model shape functions and integrate over l * + up to l * max + = 4.75 × e −1/6 GeV. This result is compatible with the NLO fit results for χ b matrix elements [10], and has been used successfully in phenomenological analyses of χ b cross sections based on the pNRQCD formalism [29].
The results for r c octet for the model shape functions are shown in fig 6. The ratio is below 1, even though dσ[χ QJ ]/dp T | NP is positive; this happens because the short-distance coefficient dσ[QQ( 3 S [8] 1 )]/dp T is smaller for the shape function formalism due to the use of a larger value of the QQ mass, which reduces the overall normalization of the shortdistance coefficient through the nonrelativistic normalization of the QQ state. The model dependence is small, amounting to less than ±0.022. Regardless of the model chosen to compute dσ[χ QJ ]/dp T | NP , the ratio r c octet is nearly flat for values of p T above about 20 GeV, consistently with the estimate we made in the beginning of this section. At lower values of p T , the ratio changes slope and steadily decreases as p T drops below from about 10 GeV. This is a combined effect of the reduction of the relative size of dσ[χ QJ ]/dp T | NP compared to dσ[χ QJ ]/dp T | octet, NRQCD , and the reduction of dσ[QQ( 3 S [8] 1 )]/dp T at m QQ = m χc compared to the m QQ = 2m case. This behavior is similar for r b octet for 1P and 2P states, as shown in fig. 7. The deviation of the ratio r b octet from one is less severe than the charmonium case. In the case of  values of p T , because the 2P state is heavier than the 1P state. Because in all cases model dependences are smaller than the usual estimates of theoretical uncertainties, we neglect the model dependence and compute dσ[χ QJ ]/dp T | NP using shape function model 1, which gives results that are close to the average of all models considered here. Model 1 seems to be a reasonable choice, since if we were to interpret the nonlocal gluonic operator vacuum expectation value E / 11 (l + ) as the probability for a gluon with momentum l to propagate to spacetime infinity, we would expect the probability to vanish as l → 0 due to confinement.
We now compute the cross sections for χ c at pp collisions at center-of-mass energy √ s = 7 TeV with rapidity cut 2 < y < 4.5. For this we need to determine the values of NRQCD matrix elements. We compute the color-singlet matrix element from the well-known result O χ c0 ( 3 P [1] 0 ) = 3Nc 2π |R (0)| 2 , which is also consistent with pNRQCD calculations of the color-singlet matrix element, and use the determination |R (0)| 2 = 0.057 GeV 5 obtained in  Figure 8. Production rate of J/ψ from χ c decays in pp collisions at √ s = 7 TeV in the rapidity range 2 < y < 4.5 computed from the shape function (red band) and NRQCD (dashed outlined band) formalisms as functions of the J/ψ transverse momentum. Experimental values (black filled circles) are obtained from the prompt J/ψ cross section and the feeddown fraction measurements from LHCb [61,65]. The "direct" experimental values (blue open squares) are obtained by subtracting the contribution from the cascade decay ψ(2S) → χ c + X followed by χ c → J/ψ + γ, which is computed from the prompt ψ(2S) cross section measurements from LHCb [66] and the measured branching fractions in ref. [67].
ref. [29] from two-photon decay rates of χ c0 and χ c2 . For the color-octet matrix element, we use the NLO best-fit value 0 ) = 0.043 at Λ = 1.5 GeV [29], as we have done for the determination of the parameters of our model shape functions. As we have stated earlier, we compute the NLO short-distance coefficients for the color-singlet and color-octet channels using the FDCHQHP Fortran package [63] with CTEQ6M parton distribution functions [58]. We set the scales for α s and the parton distributions functions to be m T = p 2 T + m 2 QQ , and vary them simultaneously between 1 2 m T and 2m T , as have been done in many phenomenological studies of quarkonium production. To match what have been measured in experiments, we compute the combined cross sections which follows from the fact that the processes χ cJ → J/ψ + γ for J = 1 and 2 are mainly E1 transitions. We take the χ cJ → J/ψ + γ branching fractions from ref. [67]. We neglect the contribution from χ c0 , because the branching fraction B χ c0 →J/ψ+γ is very small. The results for dσ[χ c → J/ψ]/dp T from NRQCD and the shape function formalisms are shown in fig. 8, compared to experimental values obtained from LHCb measurements. The experimental values for absolute cross sections have been obtained from the measurements of the ratio of χ c to J/ψ production rates in ref. [61] and the measured prompt J/ψ production rates in ref. [65]. Note that the prompt cross section measurements in ref. [65] include feeddowns from decays of ψ(2S), while our calculation does not; to facilitate a better comparison, we also show experimental values for feeddown-subtracted χ c cross sections, which we obtain from the measured ψ(2S) cross sections in ref. [66] and the branching fractions B ψ(2S)→χ cJ +X for J = 1 and 2 from ref. [67]. The difference between results from the NRQCD and the shape function formalisms are rather small and are within the theoretical uncertainties, except when p T ≈ m χc : for p T around 3-5 GeV, the cross section in the shape function formalism is smaller than the result from the NRQCD formalism, and agrees better with the feeddown-subtracted experimental values. While this difference slightly improves agreement with experiment, it only amounts to about 30%, which is the nominal size of corrections suppressed by v 2 that are neglected in both NRQCD and shape function formalisms.
We can repeat the same calculation for bottomonium: we compute cross sections for χ b (1P ) and χ b (2P ) for LHCb kinematics at √ s = 7 TeV and 2 < y < 4.5. Similarly to the case of χ c we need values for the NRQCD matrix elements. For the color-singlet matrix elements we again use the relation O χ b0 ( 3 P  Figure 10. Production rate of Υ(1S) from χ b (1P ) (left) and χ b (2P ) (right) decays times the branching fraction B ≡ B Υ(1S)→µ + µ − in pp collisions at √ s = 7 TeV in the rapidity range 2 < y < 4.5 computed from the shape function (red band) and NRQCD (dashed outlined band) formalisms as functions of the Υ transverse momentum. Experimental values (black filled circles) for p T ≤ 6 GeV are obtained by using the measured inclusive Υ cross section in ref. [68] and extrapolating the feeddown fraction measurements from LHCb [62] to smaller values of p T . The "direct" experimental values (blue open squares) are obtained by subtracting the cascade decay contributions from Υ(2S) → χ b (1P ) + X, Υ(3S) → χ b (1P ) + X, and Υ(3S) → χ b (2P ) + X, which are computed from the inclusive Υ(2S) and Υ(3S) cross section measurements from LHCb [68] and measured branching fractions in ref. [67].
calculations considered in ref. [69]. For the color-octet matrix element, we use the evolution equation in eq. (6.30) and the best-fit value for the matrix element ratio for χ c to determine O χ b0 ( 3 S [8] 1 ) at scale Λ = 4.75 GeV. On the experimental side, LHCb measurements are available for the feeddown fractions in ref. [62] for Υ(nS) from χ b (n P ) decays for n and n between 1 and 3 with n ≥ n, and the absolute cross sections for Υ(nS) as functions of p T of Υ in ref. [68]. We use the measured values for Υ(1S) cross sections and the feeddown fractions to obtain the experimental values for χ bJ (n P ) cross sections times branching fraction into Υ(1S), summed over J. To compare with the experimental values, we compute where B ≡ B Υ(1S)→µ + µ − , as functions of Υ transverse momentum p Υ T . We compute p Υ T by using We take the measured branching fractions from ref. [67]. Similarly to the charmonium case, we neglect the contribution from χ b0 on the account that its branching fraction into Υ is small. The results for dσ[χ b (n P ) → Υ(1P )]/dp T from NRQCD and shape function formalisms are shown in fig. 9, compared to experimental values. As we have done for χ c , we also show the experimental results with feeddowns from Υ(2S) and Υ(3S) subtracted by using cross section measurements in ref. [68] and branching fractions from ref. [67]. For both 1P and 2P , both formalisms yield similar results for p T 10 GeV, while the shape function formalism predicts smaller cross sections as p T decreases. For p T around 2 -4 GeV, the difference can exceed the theoretical uncertainties estimated from scale variations. The difference between the central values of the results from the two formalisms can be almost 20% for 1P and more than 30% for 2P states, which also exceed the nominal size of the order-v 2 corrections that are neglected in both formalisms. Unfortunately, the measurements for the feeddown fractions in ref. [62] are only available for p T ≥ 6 GeV, so it is currently not possible to tell which formalism would give better descriptions of smallp T behaviors of the cross sections. Since the measured feeddown fractions drop slowly as p T decreases, we could attempt to estimate the experimental values of the small-p T cross sections dσ[χ b (n P ) → Υ(1P )]/dp T at p T < 6 GeV by linearly extrapolating the feeddown fractions and multiplying by the measured Υ(1S) cross sections in ref. [68]. The result is shown in fig. 10 compared to predictions from NRQCD and shape function formalisms. It can be seen that for both 1P and 2P states, the shape function formalism leads to results that agree better with the extrapolation of the experimental results. Nevertheless, actual measurements at small p T will be needed to tell whether the shape function formalism provides better descriptions of low-p T cross sections, although this can be challenging due to diminishing experimental efficiency as p T decreases.
Before concluding this section let us discuss the possibility of further extending the formalism to even lower values of p T . In this case, there are contributions from processes where the QQ is produced with vanishing transverse momentum, and evolves into a quarkonium by emitting soft gluons. The color-singlet channel involves the contribution where QQ( 3 S is first produced at rest, which evolves into QQ( 3 P [1] J ) by emitting a gluon with momentum k. This produces a singular plus distribution in p T where the p T differential cross section diverges as p T → 0, although the integrated cross section is finite. Because the effect of gluon emission can become nonperturbative, it is tempting to replace the perturbative calculation of the gluon emission with a nonperturbative shape function, in a way similar to eq. (5.9). This way, the effect of the soft-gluon emission at small gluon momentum k would be described by the product of the QQ( 3 S [8] 1 ) cross section at vanishing transverse momentum, and the color-octet shape function with l + = 2|k|. Because the shape functions are isotropic, the shape function times the QQ( 3 S [8] 1 ) cross section predicts that the singular part of the QQ( 3 P [1] J ) cross section is isotropic in the quarkonium momentum. It is easy to see that this is not the case: the simplest way is to look into the parton cross sections for qq → QQ( 3 P J ) + g, which are available as functions of the partonic Mandelstam variables in ref. [70]. From the analytical expressions for dσ/dt, we obtain the singular parts of the differential cross sections at the parton CM frame given bŷ where we neglect the terms that are finite when the size of the quarkonium three-momentum |P | vanishes, θ is the scattering angle, and C J (cos θ) are the angular distributions given by which are normalized as +1 −1 C J (cos θ)d cos θ = 1. From the expression for the d-dimensional two-body phase space available in ref. [71], we can see that the factor |P | −1 in eq. (6.35) becomes |P | −1−2 in d = 4 − 2 dimensions, so that the cross sections contain an IR pole at |P | = 0 given bŷ Upon integrating these expressions over cos θ we reproduce the IR pole at |P | = 0 and the singular contribution at |P | → 0 given in eq. (186) of ref. [71]. Compare this with the 1 )] cross section at order α 2 s given by [71] σ[qq → QQ( 3 S which is isotropic, and reproduces the IR pole in the P -wave color-singlet cross section only after completely integrating over the scattering angle. In actual experiments, cuts are usually made so that a complete integration over the scattering angle is not possible; especially, for p T -differential cross sections p T = |P | sin θ is fixed. Hence, NRQCD factorization fails at vanishingly small values of quarkonium momentum. In contrast, at a nonzero quarkonium momentum P , there is always a finite neighborhood of P in phase space where the angles of the soft gluon momentum k can be completely integrated over, so that the cancellation of IR poles always occur. Even in this case, the size of the neighborhood of P , given by P + k, is limited by the maximum available energy |k| of the soft gluon emission. If this maximum energy is too small, such as the case where p T is smaller than the NRQCD factorization scale, radiative corrections due to soft gluon emission will involve large logarithms of the ratio of the NRQCD factorization scale and the maximum available soft gluon energy. Therefore, NRQCD factorization will become unreliable for values of the transverse momentum below the NRQCD factorization scale, and completely break down at zero transverse momentum. The shape function formalism fails for the same reason: convolving the order-α 2 s color-octet Figure 11. Tree-level Feynman diagrams for the process qq → QQ( 3 P cross section (6.37) with the perturbative shape function coincides with the singular part of the color-singlet cross section (6.35) only after completely integrating over the scattering angle, which in practice cannot be done when comparing with experiment. The root cause of the breakdown of factorization is that the infrared divergent factor in eq. (5.12) contains a dependence on the relative momentum q = (p 1 − p 2 )/2 between the heavy quark and antiquark coming from the heavy quark propagator denominator. This dependence cannot be dropped, as doing so will lead to an incorrect result for the singular contribution to the cross section in eq. (5.13). As can be seen from the Feynman diagram for the process qq → QQ( 3 P

[1]
J ) + g shown in fig. 11, the denominators of eq. (5.12) come from the virtual heavy quark lines, which have virtualities much smaller than m when the gluon momentum k is soft. Hence, this propagator must be included in the long-distance matrix element. If the QQ is produced with transverse momentum much larger than the soft gluon momentum k, we can take the average over the angles of k because the effect of variation in k on the kinematics of the quarkonium is at most order mv and can be neglected, and we obtain the usual result in eq. (5.13). However, in the case where the QQ( 3 S J ) plus a soft gluon, the gluon momentum is strongly entangled with the kinematics of the quarkonium as k = −P . In this case, the long-distance matrix element is sensitive to the details of the propagator between the heavy quark and antiquark, which does not allow a description in terms of a local operator product of the heavy quark and antiquark fields. That is, a non-local operator product of the heavy quark and antiquark fields on the same side of the projection operator is needed to describe such processes. To make matters worse, in this case it is not even clear whether it would be possible to isolate the nonlocal operator product of the heavy quark and antiquark fields into a vacuum expectation value, because soft interactions between initial and final state partons will no longer become suppressed, unlike the case of large-p T production. A similar conclusion was obtained in the analysis of transverse-momentum dependent factorization of P -wave quarkonium production in ref. [72], where the authors considered the process gg → QQ( 3 P [1] J ) + X, which is more complicated because the process gg → QQ( 3 S [8] 1 ) vanishes at tree level.
As we have shown, this breakdown of NRQCD factorization at p T 0 signals a failure of the factorization formalism in terms of local operator matrix elements. This reveals an important distinction between the description of B hadron production in heavy quark effective theory [73] and quarkonium production in NRQCD. In heavy quark effective theory, where long-distance matrix elements involve the b quark field in place of the QQ bilinears, the b quark field is always a local field. In contrast, in quarkonium production nothing prevents the heavy quark and antiquark from being produced nonlocally; this nonlocality can be neglected if it is disentangled from the short-distance process, and in this case we recover the standard form of NRQCD factorization. However this is not always guaranteed to happen, and can fail for production with small transverse momentum as we have seen from the example above.

Summary and conclusion
In this work we investigated the application of the shape function formalism for resumming kinematical effects in inclusive production rates of P -wave heavy quarkonia. In this formalism, the products of short-distance coefficients and NRQCD matrix elements are replaced by convolutions of short-distance coefficients and shape functions. The shape functions encode the kinematical effects arising from the motion of the heavy quark and antiquark pair in the heavy quarkonium rest frame. They are nonperturbative functions of a lightlike momentum l, which corresponds to the difference between the momentum of the quarkonium and the momentum of the heavy quark-antiquark pair. While the formalism was originally developed in refs. [21,22] based on tree-level calculations of quarkonium production, we have extended the application of the formalism to next-to-leading order in the strong coupling, which requires treatment of the ultraviolet divergences that arise from renormalization of NRQCD matrix elements. Due to the mixing between color-singlet and color-octet channels induced by the infrared divergence in the color-singlet P -wave production cross section, the short-distance coefficients differ between the NRQCD and shape function formalisms. Based on the perturbative calculations of the shape functions at next-to-leading order accuracy in sec. 3, we have derived in sec. 5 the matching conditions that are valid at one-loop level which allow determination of the short-distance coefficients in the shape function formalism from the standard NRQCD ones. This allows calculation of the cross sections in the shape function formalism at one-loop level using the known NRQCD short-distance coefficients in the MS scheme, once the nonperturbative shape functions are known.
In sec. 4 we presented an improved analysis of the nonperturbative shape functions for P -wave heavy quarkonia. By using the potential NRQCD formalism developed in refs. [14,16,28,29], we reproduced the proposed form of the color-singlet shape function in ref. [21] that is proportional to δ(l + ), where l + is the + component of the lightlike momentum l associated with the shape function formalism. Together with the requirement that the normalization of the shape function must be equal to the corresponding NRQCD matrix element, this completely fixes the color-singlet shape function up to corrections suppressed by powers of v. In the case of the color-octet shape function, in existing studies it was assumed that it would not be possible to calculate them except within models. While first-principles determinations of the nonperturbative color-octet shape functions may be currently out of reach, we find that the large-l + asymptotic form fixed by the perturbative calculation, and the normalization condition that its integral over l + must reproduce the NRQCD matrix element, strongly constrain the l + dependence of the shape function. Furthermore, we showed through a nonperturbative potential NRQCD calculation of the shape functions that the l + dependence of the color-octet shape functions must be same for all P -wave heavy quarkonia, independently of the radial excitation or the heavy quark flavor, up to corrections suppressed by powers of v. This strongly constrains the behavior of the shape function, and the model dependence is constrained in the small l + region. As the contribution from the small l + region to the cross section is suppressed, this model dependence has only a small effect on the P -wave quarkonium cross section. Hence, the analysis in this work diminishes the model dependence of the shape function formalism for production of P -wave heavy quarkonia.
Based on the results for the one-loop perturbative matching conditions and nonperturbative shape functions obtained in this work, we presented two phenomenological applications. First, in section 6.1 we used the fact that the color-octet shape function describes soft-gluon emissions to match the accuracies in α s of the singular radiative corrections near threshold between the color-singlet and color-octet states. This alleviates the problem in fixed-order calculations that, when the short-distance coefficients are computed to same accuracy in the strong coupling, the singular threshold effects are truncated to different orders in α s in the color-octet and color-singlet channels. This can be problematic especially at very large transverse momentum, where the threshold effects can be enhanced. As we have shown in section 6.1, standard NRQCD calculations can yield negative production rates of χ c ; the situation is much more improved when we match the singular threshold effects between the two channels using the results from the shape function formalism, and the positivity of the cross sections are ensured to very large values of transverse momentum. While this does not eliminate the need for threshold resummation, the matching of the threshold effects considered in this work can be useful for consistent treatment of threshold resummation throughout the singlet and octet channels.
In section 6.2, we investigated the impact of the nonperturbative shape function on the transverse-momentum dependent cross sections of χ c and χ b . We employed a variety of model functions for the color-octet shape function to compute the nonperturbative corrections, where the model parameters are constrained by the phenomenologically obtained values of the color-octet matrix elements. As we have argued above, the model dependence only affects the l + → 0 behaviors of the shape functions, and have little effect on the cross sections, so that the model dependence can even be neglected compared to other theoretical uncertainties. We find that at large p T , the effect of nonperturbative corrections is mainly to change the overall normalization of the color-octet contribution, which is inconsequential in phenomenology because a change in normalization can be reabsorbed into the color-octet matrix element. This situation gradually changes as we go to lower transverse momentum, and at values of the quarkonium transverse momentum similar or smaller than the heavy quarkonium mass, the cross section in the shape function formalism undershoots the NRQCD prediction by more than the usual size of estimated theoretical uncertainties obtained from scale variations. Although this result is interesting, it is currently not possible to tell whether the shape function formalism gives a better description of the P -wave heavy quarkonium cross sections at small p T ; in the case of χ c , the difference between the two formalisms is not too large compared to the nominal size of relativistic corrections suppressed by v 2 , which is usually taken to be about 0.3 for charmonium. In the case of χ b , the results from the shape function formalism does tend to agree better with experimental values obtained by extrapolating the feeddown fractions of Υ from χ b decays to p T ≤ 6 GeV; nevertheless, actual measurements will be necessary to tell if this is indeed the case.
As we have mentioned, the P -wave quarkonium production cross sections are also important for understanding the feeddown contributions in S-wave quarkonium production rates. It would also be interesting to apply the shape function formalism for production of S-wave quarkonia, such as J/ψ, ψ(2S), and Υ states. Production of S-wave quarkonia is phenomenologically more significant, and measurements are available for a wider range of transverse momentum. As have been discussed in ref. [16], the color-octet NRQCD matrix elements that appear in S-wave quarkonium production rates not only involve logarithmic UV divergences that induce mixing between different color-octet channels, but also quadratic power divergences as well, unlike the P -wave case. We can expect that these quadratic power divergences will appear in color-octet shape functions as contributions proportional to l + . Since in perturbative QCD these contributions will lead to scaleless power UV divergences, they will need to be subtracted in order to make contact with NRQCD short-distance coefficients that are computed in dimensional regularization. Due to this the nonperturbative shape functions may affect the S-wave cross sections differently from P -wave cross sections. It would be interesting to see how the application of the shape function formalism will affect S-wave quarkonium cross sections.