QCD corrections to J/psi polarisation in pp collisions at RHIC

We update the study of the polarisation of J/psi produced in proton-proton collisions at RHIC at sqrt(s)=200 GeV using the QCD-based Colour-Singlet Model (CSM), including next-to-leading order partonic matrix elements from gluon and light quark fusion and leading-order contributions from charm-quark initiated processes. To do so, we also evaluate the corresponding cross section differential in P_T which agrees qualitatively with the measurements of PHENIX in the central and forward regions at low P_T -- for instance below 2 GeV --, while emphasising the need for Initial State Radiation (ISR) resummation. At mid P_T, we also compare the measurements from PHENIX and STAR with the same evaluation complemented with the dominant alphaS^5 contributions (NNLO*). We find a reasonable agreement with the data. Regarding the polarisation, as shown for previous studies at larger sqrt(s) and P_T, the polarisation pattern from gluon and light quark fusion in the helicity frame is drastically modified at NLO and is shown to be increasingly longitudinal. The yield from charm-gluon fusion is found to be slightly transversally polarised. Combining both these contributions with a data-driven range for the polarisation of J/psi from chi_c, we eventually provide an evaluation of the polarisation of the prompt J/psi yield which is in a good agreement with the experimental data from PHENIX both in the central and forward regions.


Introduction
Until recently, the numerous puzzles in the prediction of quarkonium-production rates at hadron colliders were attributed to non-perturbative effects associated with channels in which the heavy quark and antiquark are produced in a colour-octet state [1,2,3,4]. It is now widely accepted that α 4 s and α 5 s corrections to the CSM [5] are fundamental for understanding the P T spectrum of J/ψ and Υ produced in high-energy hadron collisions [6,7,8,9,10,11,12]. The effect of QCD corrections is also manifest in the polarisation predictions. While the J/ψ and Υ (commonly denoted Q hereafter) produced inclusively or in association with a photon are predicted to be transversally polarised at LO, it has been recently emphasised that their polarisation at NLO is increasingly longitudinal when P T gets larger [8,10,13,14].
In a recent work [15], we have also shown that hard subprocesses based on colour singlet QQ configurations alone are sufficient to account for the observed magnitude of the P T -integrated cross section. In particular, the predictions at LO [5] (Fig. 1 (a)) and NLO [6,7,8] (Fig. 1 (b,c,d)) accuracy are both compatible with the measurements by the PHENIX collaboration at RHIC [16,17] within the present uncertainties. This provided some indications that the computations are carried in a proper perturbative regime. This agreement is improved when hard subprocesses involving the charm-quark distribution of the colliding protons are † Permanent address taken into consideration. These constitute part of the LO (α 3 S ) rate ( Fig. 1 (e)) and are responsible for a significant fraction of the observed yield, as we have argued in [15].
In this Letter, we proceed to the evaluation of the P T dependence of the polarisation of the J/ψ produced at RHIC, both in the central and forward rapidity regions. In the section 2, we expose the procedure used to compute the yield and the polarisation at NLO (up to α 4 S ) from gg & gq fusion and from cg fusion at LO (at α 3 S ) as well as the procedure to obtain a first evaluation of the leading P T contributions from gg and gq fusion at α 5 S (NNLO ). In the section 3, we present our results. First, we show the yields differential in the rapidity y and P T along with the polarisation vs P T from gg & gq and from cg fusion separately. Then we show our results when they are combined. Finally we combine the latter predictions for the direct yield polarisation with an essentially data-driven estimation of the polarisation for J/ψ from χ c and we compare our results with the PHENIX data in both rapidity ranges. The last section gathers our conclusions.

Cross section
In the CSM [5], the matrix element to create a 3 S 1 quarkonium Q of momentum P and polarisation λ accompanied by other partons, noted j, is the product of the amplitude to create the corresponding heavy-quark pair, a spin projector N(λ|s 1 , s 2 ) and R(0), the radial wave function at the origin in the configuration space, obtained from the leptonic width [18], namely where P = p Q + pQ, p = (p Q − pQ)/2, s 1 ,s 2 are the heavyquark spin and δ ii / √ N c is the projector onto a coloursinglet state. In the non-relativistic limit, N(λ|s 1 , s 2 ) can be written as 2m Qv ( P 2 , s 2 )γ µ u( P 2 , s 1 ) where ε λ µ is the polarisation vector of the quarkonium. The sum over the quark spin yields to traces evaluated in a standard way.
In our evaluation, we use the partonic matrix elements from Campbell, Maltoni and Tramontano [6] to compute the LO and NLO cross sections from gluon-gluon and lightquark gluon fusion. We guide the reader to [6] for details concerning the derivation of the cross section at α 4 S , the corresponding expressions at α 3 S can be found in [19]. In the case of the cg fusion (at LO), we use the framework described in [20] based on the tree-level matrix element generator MADONIA [21].
Moreover, in order to illustrate the expected impact of NNLO QCD corrections for increasing P T , we also present the results when the leading-P T contributions at α 5 S evaluated along the lines of [10] are added to the previous contributions. At this order, the last kinematically-enhanced topologies open up. These exhibit a P −4 T fall off of their differential cross section in P 2 T , typical of a single particle exchange in the t channel. Going further in the α S expansion cannot bring any further kinematical enhancement as regards the P T dependence. As a consequence, above α 5 S , usual expectations for the impact of QCD corrections would then hold and these should be then well taken into account by a K factor multiplying the yield at NNLO accuracy, which would be independent of P T and of a similar magnitude as those of other QCD processes. In other words, a further cross-section modification between the NNLO and N 3 LO results by an order of magnitude would not be acceptable at any P T , while it is expected to be so between the cross sections at LO (P −8 T ), NLO (P −6 T ) and NNLO (P −4 T ) for large enough P T , simply owing to their different P T scalings.
The procedure used here to evaluate the leading-P T NNLO contributions is exactly the same as in [10]. Namely, the real-emission contributions at α 5 S are evaluated using MADONIA by imposing a lower bound on the invariantmass of any light partons (s i j ). For the new channels opening up at α 5 S with a leading-P T behaviour, and which specifically interest us, the dependence on this cut is to get smaller for large P T since no collinear or soft divergences can appear there. For other channels, whose Born contribution is at α 3 S or α 4 S , the cut would produce logarithms of s i j /s min i j , which are not necessarily small. Nevertheless, they can be factorised over their corresponding Born contribution, which scales as P −8 T or P −6 T , and are thence suppressed by at least two powers of P T with respect of the leading-P T contributions (P −4 T ). The sensitivity on s min i j is thus expected to come to nothing at large P T . This argument has been checked at α 4 S (NLO vs. NLO ) for Υ [10] and ψ [12] in the inclusive case as well as in association with a photon [14].

Polarisation
The polarisation parameter α (also called λ) is computed by analysing the distribution of the polar angle θ between the + direction in the quarkonium rest frame and the quarkonium direction in the laboratory frame. This definition of θ is referred to the analysis of the polarisation (or spin-alignement) in the helicity frame 2 .
The normalised angular distribution I(cos θ) reads from which we can extract α bin by bin in y or P T . The interpretation of α is immediate when one relates it to the polarised cross sections: For σ T σ L , thus for a yield purely transversally polarised, α 1 while for σ L σ T (a yield purely longitudinally polarised), α −1. On the contrary, if no direction is favoured, one expects σ T = (σ T x + σ T y ) = 2σ L , which corresponds to α 0, and thus no θ dependence. We emphasise here that the relations Eq. (2) and Eq. (3) do not depend on the definition chosen for θ (the frame definition). Yet, the results obtained do depend on it: a yield transversally polarised in one frame can be longitudinal in another. In particular, the expressions that one would obtain for σ L and σ T change from one frame to another.

gg and gq channels
We first present the results from the gg and gq channels at LO (α 3 S ) and NLO (α 3 S +α 4 S ) in terms of differential cross sections as function of y and P T , which we compare to the experimental data for prompt J/ψ multiplied by the expected fraction of direct J/ψ (59 ± 10)% (see [15]).
As discussed in [15], the experimental data for the P T integrated cross section from PHENIX [16,17] are qualitatively well reproduced by the highest value of the theoretical bands, both at LO and NLO (see Fig. 2). Let us also note that choosing a lower value than 1.4 GeV for the charm quark mass -along the same lines as what has been done in studies of charm production at RHIC [25]-would lead to higher cross sections. However a part of the increase of the cross section could be attributed to the approximation M J/ψ = 2m c , which would then induce an artificially low bound-state mass in the kinematics. Such an effect is not at LO and NLO accuracy compared to the PHENIX data [16,17]. The theoretical-error bands for LO and NLO come from combining the uncertainties resulting from the choice of µ f , µ r , m q . We have used the LO set CTEQ6 L [24] for the yield at LO, the NLO set CTEQ6 M [24] for the yield at NLO.
easily quantifiable without a dedicated study. Therefore, we do not use values lower than 1.4 GeV for m c here. Before discussing the results of the cross section differential in P T at LO and NLO, it is important to note that we are working in the low P T region where the perturbative expansion is not always reliable. While the comparison between the results at LO and NLO integrated in P T shown on Fig. 2 gives some good indication that the perturbative expansion works well, it may not be so when one considers the cross section differential in P T . First, the yield at NLO has shown negative -thence unphysical-values in some bins in P T and y. When this occurs, this is always for the lowest P T bin. This shows that the virtual corrections at α 4 S (with a negative relative sign) are large. This may also be the case for the virtual corrections at α 5 S which are currently unknown and possibly with an opposite sign to the ones at α 4 S . This is a known issue which can be partly solved by resumming Initial State Radiation (ISR). Such a resummation has been carried out for the Υ production at the Tevatron using the Colour Evaporation Model [26]. One of the main outcomes of the latter study is that the P T dependence of the cross section is significantly affected at low P T (up to roughly m Υ /2) by ISR. A similar impact of ISR is expected for the charmonium family, albeit for a narrower range in P T .
Since we have not carried out this resummation 3 , the results that we have obtained for P T below m J/ψ /2 for dσ/dP T shown on Fig. 3 (as well as later for the polarisation) are to be taken with a grain of salt. That being said, one sees on Fig. 3 that the slope for the yield at NLO is milder than for the LO, though too steep to reproduce the data from PHENIX [16] and STAR [27]. Yet, the general effect of ISR is to lower the cross section at very low P T and to rise it a little up to P T 1 − 2m c .
In addition, one expects a significant contribution from α 5 S contributions at larger P T . A complete evaluation of NNLO corrections is not yet available. For now, we can only rely on a study of the leading-P T NNLO contributions (NNLO ) as done in [10,12]. However, for such low values of P T , the latter method is a priori not reliable and we could not extend its evaluation (red band) in Fig. 2 (a) 4 below P T = 5 GeV. Overall, the comparison presented allows us to think that a more complete analysis (through ISR resummation and matching with leading P T contributions for instance) could show that the P T dependence given by the CSM agrees with the experimental measurements in the low and mid P T regions at RHIC energies. We now turn to the discussion of the polarisation results vs P T (Fig. 4) in the helicity frame. As can be seen on Fig. 4, the complete modification of the polarisation pattern between the LO and the NLO results observed in previous works [8,10,13,14] is confirmed, despite the error band  In the forward region, the yield at NLO is transverse at low P T and becomes longitudinal as soon as P T ≥ m c .
As discussed previously, the result obtained in the lowest P T bin (below 1 GeV) are likely to be subject to higher QCD corrections and the band should not be extrapolated down to P T = 0. Second, the extraction of the yield polarisation at NLO at low P T is highly demanding in terms of computer time. Typically, the evaluation of the yield polarisation in a single P T bin of 500 MeV with statistical fluctuation less than 20 % requires O(10 8 ) numerical evaluations of the integrand. These fluctuations add up to the usual theoretical uncertainties in the bands shown on Fig. 4.

cg channel
We first present the results from the cg LO (α 3 S ) in terms of differential cross sections as function of y and P T , which we compare to the experimental data for prompt J/ψ multiplied by the expected fraction of direct J/ψ.
In [15], we have computed the yield integrated over P T from cg channels and used the LO set CTEQ6.5c [28] based on a recent global PDF fit including Intrinsic Charm (IC). More precisely, we have compared the results for three choices for the charm distribution: (i) without IC [c(x, µ 0 ) = 0 (µ 0 = 1.   ( x c+c = 2.4%). In the following, we shall carry on our study of the polarisation with the sealike IC. The result are qualitatively the same with other choice for c(x). By comparing the data to our results for the contribution from cg fusion shown on Fig. 5, one sees that it accounts for 5 to 40 % of the observed yield depending on the usual theoretical uncertainties. Comparable fractions are obtained for other c(x) as discussed in [15].
Regarding the P T dependence shown in Fig. 6 for both rapidity regions, one observes that the contribution is falling too fast compared to the PHENIX data. Nevertheless, the appearance of NLO contribution such as cg → J/ψcg (e.g. Fig. 1 (f)) at larger 5 P T , and to a lesser extent the effect of ISR at small P T , is expected to smear the curve out to larger P T . 5 These were for instance analysed in the fragmentation approximation in [30].
To what concerns the polarisation pattern (Fig. 7), in the central region, it is similar to the one computed for J/ψ + cc in [7], namely nearly unpolarised as soon as P T ≥ m c , although not completely unpolarised with α 0.2 for mid P T . In the forward region, it becomes somewhat more transversally polarised and shows a dependence on the mass and scale choices as indicated by the widening band.

Polarisation of the direct yield
In order to obtain the polarisation of the direct yield, one has to combine the polarisation from both contributions taking into account their proper weight. Schematically one has: where ∆σ is the differential cross section for one given process integrated in a given bin in y and P T . 5 2) regions at NLO + (and cg LO + NNLO for (a)) compared to the PHENIX [16] and STAR [27] data. The theoretical-error bands come from combining the uncertainties resulting from the choice of µ f , µ r , m q , see text. : Comparison between α(P T ) for J/ψ directly produced at NLO + (see text) and the PHENIX data [31,32] in both rapidity regions in pp at √ s NN = 200 GeV.
Beforehand, we present results for the NLO + yield, namely the yield at NLO accuracy from gg and gq fusion added to the yield from cg fusion at LO accuracy 6 . The sum of both contributions differential in P T is compared to the PHENIX and STAR data on Fig. 8 a) and b). In the central region, the yield at NNLO from gg and gq fusion is also shown (with cg at LO added). The computation being close to the data for P T < 4 GeV, it is reasonable to compare them to the PHENIX measurements [31,32] as done on Fig. 9.
Except for the lowest P T point in the central region, the direct NLO + yield polarisation is compatible with the 6 The NLO corrections to cg → J/ψX are not yet known.
(prompt) data 7 , indicating a small impact of the ψ(2S ) and χ c feed-downs on the polarisation. The same conclusion holds from the s-channel cut analysis [33,34], where a good agreement with experimental data was also obtained, at least in the central region. It is worth recalling that the latter analysis was done by neglecting the usual contribution from the CSM which are the purpose of this work and are evidently not negligible. Since the results for the polarisation are similar in both analysis, a combined study would follow the same trend.
Finally, the disagreement between the data and the present evaluation in the lowest P T bin is not worrisome in view of the expected impact of higher QCD corrections in this bin.

Data-driven evaluation of the feed-down effect on the polarisation
For the time being, there exists no measurement of the direct J/ψ polarisation after extraction of both the B and χ c feed-downs. While at RHIC energy the former contributes less than 5 % of the full yield [35] and thence should not affect its polarisation, the χ c feed-down (up to 30-40 %) may impact strongly on the observed values of α. It is however possible to constrain its effects by using existing data from Hera-B on σ χ c1 /σ χ c2 for instance, relying on the dominance of the E1 transition between the χ c 's and the J/ψ.
Indeed, using E1 dominance [36], one can obtain a first approximation of the yield of longitudinally (transversally) polarised J/ψ, denoted N h=0 ) from the parent χ c , in terms of simple 8 relations involving the polarised χ c yields, these are Br( 3 P 2 → 3 S 1 γ) 8 We should emphasise that these relations neglect effects that could arise from the difference in the definition of the helicity frame of the χ c and that of its decay product (the J/ψ here). These effects are expected to vanish for large momenta for the quarkonium. Such an approximation may not be accurate for the smallest P T 's considered here. 200 GeV using a sealike charm distribution compared to the PHENIX data [16]. The theoretical-error bands for LO and NLO come from combining the uncertainties resulting from the choice of µ f , µ r , m q . and N |h|=1 Let us analyse now the two otherwise extreme cases α max from χ c and α min from χ c : • α max from χ c is reached when N |h|=1 is maximised and N h=0 minimised. This happens when the 3 P 1 yield is fully h = 0 and the 3 P 2 one is fully |h| = 2. This indeed gives N h=0 3 S 1 = 0 and thus α max from χ = +1.
• α min from χ c is reached when N |h|=1 is minimised and N h=0 maximised. This happens for a 3 P 1 yield fully |h| = 1 and the 3 P 2 one fully h = 0. This gives: Recently, the Hera-B [37] and CDF [38] collaborations have measured the ratio of R 12 = σ χ c1 Br(χ c1 →J/ψγ) σ χ c2 Br(χ c2 →J/ψγ) which we can therefore use in Eq. (5). It is however worth noticing that both results are somewhat different, R 12 = 1.0 ± 0.4 for Hera-B for 0 < P T < 2 GeV at √ s = 41.6 GeV and R 12 = 2.5 ± 0.1 for CDF for P T > 4.0 GeV at √ s = 1.96 TeV. Since we focus here on low P T data, we prefer to opt for the Hera-B value, which gives for the central value Supposing that about 30% of the J/ψ come from χ c decays independent of P T in the range considered here, we expect a partial contribution to the polarisation ranging from 0.3 × (+1) to 0.3 × (−0.47). Regarding the other 70%, one can simply multiply the result obtained above for the direct yields (section 3.3) by 0.7, since in a good approximation the polarisation of J/ψ from ψ(2S ) is expected to be identical to the direct one.

Comparison with the PHENIX data
Combining the result obtained for the direct yield at NLO + and the extreme expected polarisations of J/ψ from χ c , one can compare (Fig. 10) our results for the CSM prompt J/ψ yield with the PHENIX data in the central [31] and the forward [32] regions. Within the present admittedly large theoretical uncertainties, the agreement with the data is good. 7 α for the direct NLO + (two dashed lines) and the PHENIX measurements in both rapidity regions [31,32].
Three comments are in order: a) the result should not be straightforwardly extrapolated to P T = 0, b) at larger P T , α 5 S (for gg + gq) and α 4 S (for cg) will eventually dominate due to their P −4 T scaling. While the former is expected from similar studies [10,11,12] to be longitudinal, we have qualitatively checked that the yield from the latter (more precisely from cg fusion at NLO ) is similar to cg → J/ψc at LO; c) for BHPS c(x) (showing a peak of the charm distribution around x ≈ 0.1 − 0.3), one expects a relative enhancement of the yield from cg → J/ψX for increasing P T , i.e. for a momentum fraction of the charm quark in the proton approaching the value where c BHPS (x) peaks. Yet, a dedicated study of these effects is needed to draw more quantitative conclusions. Indeed, depending on the relative yield of gg & gq fusion at NNLO and cg fusion at LO (or NLO ) at large P T , the polarisation may end up to be strongly longitudinal or slightly transversal. In any case, as opposed to nonrelativistic QCD analyses of including colour octets [39] (see [40] for a recent application for RHIC energies) transitions, the CSM yield is not expected to become strongly transverse for increasing P T .
Our results for gg and qg fusion are in qualitative agreement with the analysis carried out by Baranov and Szczurek [41] using the k t factorisation approach. However the latter study predicts a dominance of χ c2 feed down over the prompt J/ψ yield which may conflict both with the preliminary upper value [35] for the χ c feed down (< 42% at 90% C.L.) by the PHENIX collaboration and with the Hera-B results for the ratio σ χ c1 /σ χ c2 . We further note that some contributions appearing at α 5 S and enhanced by log(s) [42] also produce a yield which is mainly longitudinal.

Conclusions and outlooks
In conclusion, we have carried out the first NLO analysis in the Colour-Singlet Model of the polarisation of the direct J/ψ yield at RHIC including contributions from c-quarkgluon fusion. Our result for the yield differential in P T is in near agreement with the measurement at low and mid P T both in the central and forward rapidity regions. As regards the polarisation, our evaluation for the direct yield is in good agreement with the PHENIX data for prompt J/ψ. This therefore points at a small impact of the feed-downs on this observable in this kinematical region.
Using constraints from existing data, namely the fraction of χ c feed-down and the ratio σ χ c1 /σ χ c2 , we have determined likely extreme values of the polarisation of the J/ψ from χ c by relying on the E1 dominance of the transition χ c → J/ψγ. This has enabled us to extrapolate our polarisation evaluation of the direct yield to the prompt one. The results obtained are also in good agreement with the PHENIX data.
Motivated by this agreement, the next step would be to study the Cold Nuclear Matter (CNM) effects on the yield polarisation in the CSM framework. These effects are expected to come mainly from the shadowing of parton distributions and the final-state interactions between the cc pair and the CNM. We could for instance take benefit of the Glauber framework JIN [43] used in [44] to study the CNM effects on the J/ψ yield in dA and AA collisions at RHIC; the latter can indeed deal with the polarisation information. Such a study would surely be very expediently done in order to refine the arguments presented in [45] to employ quarkonium polarisation in heavy-ion collisions as a possible signature of the quark-gluon plasma.