Transverse Single Spin Asymmetry in p+p^\uparrow \rightarrow J/ψ+X
Rohini M. Godbole, Abhiram Kaushik, Anuradha Misra, Vaibhav Rawoot, Bipin Sonawane
CCERN-TH-2017-047
Transverse Single Spin Asymmetry in p + p ↑ → J/ψ + X Rohini M. Godbole
Centre for High Energy Physics, Indian Institute of Science, Bangalore, India. ∗ Abhiram Kaushik
Centre for High Energy Physics, Indian Institute of Science, Bangalore, India. † Anuradha Misra
Department of Physics, University of Mumbai, Mumbai, India ‡ Vaibhav Rawoot
Department of Physics, University of Mumbai, Mumbai, India. § Bipin Sonawane
Department of Physics, University of Mumbai, Mumbai, India. ¶ (Dated: June 28, 2018) Abstract
We present estimates of transverse single spin asymmetry (TSSA) in p + p ↑ → J/ψ + X withinthe colour evaporation model of charmonium production in a generalized parton model (GPM)framework, using the recently obtained best fit parameters for the gluon Sivers function (GSF)extracted from PHENIX data on TSSA in p + p ↑ → π + X at midrapidity. We calculate asymmetryat √ s = 200 GeV, and compare the results with PHENIX data on TSSA in the process p + p ↑ → J/ψ + X . We also present estimates for asymmetry at √ s = 115 GeV corresponding to the proposedfixed target experiment AFTER@LHC and at √ s = 500 GeV corresponding to the higher RHICenergy. Finally, we investigate the effect of the transverse momentum dependent (TMD) evolutionof the densities involved, on the asymmetry. ∗ Electronic address: [email protected] † Electronic address: [email protected] ‡ Electronic address: [email protected] § Electronic address: [email protected] ¶ Electronic address: [email protected] a r X i v : . [ h e p - ph ] D ec . INTRODUCTION The role of intrinsic transverse momentum and spin distribution of partons inside thehadrons in explaining the azimuthal asymmetries arising in experiments involving polarizedbeams [1–3] has been a subject of keen interest during the past decade. Transverse SingleSpin Asymmetries (TSSA’s), measured in meson production in hadronic collisions as well asin semi inclusive deep inelastic scattering (SIDIS) experiments, provide useful informationtowards understanding the transverse dynamics of partons inside the hadrons.It is known for long [4] that the conventional QCD factorization with collinear partondistribution functions (PDFs) and fragmentation functions (FFs) is not sufficient to accountfor the single spin asymmetries (SSAs) observed experimentally. It is now well establishedthat the SSAs observed in the hadroproduction of mesons and semi inclusive deep inelasticscattering (SIDIS) [5–8] can be explained in terms of orbital motion of quarks and gluonsand the spin structure of nucleons. One of the two approaches towards building up asuitable framework involves generalization of the concept of collinear PDFs by including thetransverse momentum and spin dependence of the partons into PDFs and FFs, which arethen collectively referred to as TMDPDFs. This approach known as TMD approach has beensuccessful in explaining the existing data in some processes, for example in inclusive pionproduction in pp collisions [9] and Z -boson production in pp and p ¯ p collisions [10]. In theTransverse Momentum Dependent factorization scheme [11–14], the hadronic cross-sectionis then expressed as a convolution of TMDPDFs, TMD FFs and partonic cross-section.One of the interesting TMDPDFs, which has been a subject of a large number of stud-ies is the Sivers Function, which can be interpreted as the number density of unpolarizedquarks or gluons inside a transversely polarized proton. Information on the Sivers functioncan shed light on the 3-dimensional structure of nucleon and may also provide an esti-mate of the orbital angular momentum of quarks and gluons [15, 16]. It was introducedby Sivers in Ref. [17, 18] wherein a transverse momentum dependent PDF, now knownas Sivers function, was introduced for the first time to account for the large asymmetriesobserved in pp ↑ scattering [17]. Quark Sivers function, which is usually coupled to rea-sonably well parametrized unpolarized fragmentation functions, has been one of the firstTMDs extracted from data [15]. All the known information on quark Sivers function hasbeen obtained from SIDIS data and γ ∗ − N scattering in the centre of mass frame. Initial2xtractions of up and down quark Sivers functions, from HERMES and COMPASS data,were obtained using parameterizations that did not take into account sea quarks [19, 20].Later, magnitude of asymmetry for K + measured at HERMES indicated the possible contri-bution of sea quarks and fits were obtained taking into account contributions of sea quarksalso [19]. Parametrizations of Sivers function for up and down quarks have been foundto be in agreement with light-cone models [21] and a quark-diquark spectator model [22].Extractions of quark Sivers function taking into account TMD evolution have also beenobtained by performing a global fitting of all data on Sivers asymmetry in SIDIS from HER-MES, COMPASS and Jefferson Lab [10] and based on these, predictions have been made forSivers asymmetry in the DY process and W -boson production. Although there are manyparameterizations available for quark Sivers function, not much information is available ongluon Sivers function (GSF). Some of processes which have been proposed for obtaining in-formation about GSF are pp ↑ → γ + X [23], pp ↑ → D + X [24, 25] , pp ↑ → γ + jet + X [23, 26], pp ↑ → γ ∗ + X → µ + µ − + X [23] and pp ↑ → η c/b + X [27].Heavy flavour production- both open and closed- are considered to be clean probes ofthe GSF since heavy quarks are predominantly produced via gluon-gluon fusion and thuscan be used to isolate gluon dynamics within hadrons [24, 28–32]. The possibility of get-ting information on the GSF by looking at D -meson production in polarized proton-protonscattering at RHIC has been discussed by Anselmino et al. [24] using a saturated GSF.The first phenomenological study on the gluon Sivers function has recently been performedby D’Alesio, Murgia and Pisano [33], wherein the gluon Sivers function has been fitted tomidrapidity data on transverse SSA in pp → π + X measured by PHENIX collaboration atRHIC [1]. In our previous work, where we made predictions for TSSA in electroproductionof J/ψ assuming a transverse momentum dependent factorization within colour evaporationmodel (CEM) of charmonium production, we had used parameterization of the gluon Siversfunction suggested by Boer and Vogelsang [34], in which the x -dependence of the GSF wasmodeled on that of the u and d quark Sivers functions. We will call these BV parametersin the following. In our recent work on TSSA in D -meson production [35], we have usedthe directly fitted GSF parameters of Ref. [33] and have compared the estimates using thesewith the results obtained using BV parameters.In this work, we present predictions for TSSA in the process pp ↑ → J/ψ + X using recentdirectly fitted parameters, which we will call DMP fits [33]. We compare these results with3redictions obtained using the BV parameters, which are based on experimentally fittedquark Sivers parameters [34]. We then compare our results with the recent measurementsof Sivers asymmetry at PHENIX experiment in J/ψ production [36]. A similar comparisonwith PHENIX results has been performed recently in Ref. [32] using colour singlet model(CSM) and a maximized GSF, whereas in our work we use CEM and compare variousparameterizations of GSF available. We also present estimates of asymmetry for futureproposed experiments at AFTER@LHC which is a fixed target experiment with √ s =115 GeV and for √ s = 500 GeV which will be explored at RHIC. Finally, to assess theeffect of QCD evolution on asymmetry, we compare the predictions based on DGLAP andTMD evolution of the unpolarized TMDPDF and gluon Sivers function. This comparisonis performed using the BV parameters because direct fits of the GSF with TMD evolutiontaken into account are currently not available. II. FORMALISMA. TSSA in the process p + p ↑ → J/ψ + X using colour evaporation model We consider the transverse single spin asymmetry, A N = dσ ↑ − dσ ↓ dσ ↑ + dσ ↓ (1)for the process p + p ↑ → J/ψ + X using the colour evaporation model [37] of J/ψ production.In the colour evaporation model, the total cross-section for the production of
J/ψ atleading order (LO) is proportional to the rate of c ¯ c production integrated over the invariantmass-squared range 4 m c to 4 m D , where m c is the mass of the charm quark and m D is theopen charm threshold [38]: σ p + p → J/ψ + X = F J/ψ (cid:90) m D m c dM c ¯ c (cid:90) dx a dx b (cid:20) f g/p ( x a ) f g/p ( x b ) d ˆ σ gg → c ¯ c dM c ¯ c + (cid:88) q f q/p ( x a ) f ¯ q/p ( x b ) d ˆ σ q ¯ q → c ¯ c dM c ¯ c (cid:21) , (2)where q = u, d, s, ¯ u, ¯ d, ¯ s .Here, the CEM parameter F J/ψ is the fraction that gives the probability of
J/ψ productionbelow D ¯ D threshold. 4ere, we use a phenomenological approach referred to in literature as the GeneralizedParton Model (GPM), which has been used to estimate SSAs in several processes like pp ↑ → D + X [24, 32, 35], pp ↑ → γ + jet + X [23, 26] and pp ↑ → π + X [9] for whichTMD factorization has not yet been established. A rigorous treatment will require inclu-sion of intrinsic transverse momentum effects through a consideration of higher twist effects.However, motivated by the phenomenological successes of the GPM [9, 39], we assume ageneralization of CEM expression and include TMDPDFs thus expressing the cross sectionfor the transversely polarized scattering process p + p ↑ → J/ψ + X as σ p + p ↑ → J/ψ + X = F J/ψ (cid:90) m D m c dM c ¯ c (cid:90) dx a d k ⊥ a dx b d k ⊥ b (cid:40) f g/p ( x a , k ⊥ a ) f g/p ↑ ( x b , k ⊥ b ) d ˆ σ gg → c ¯ c dM c ¯ c + (cid:88) q (cid:20) f q/p ( x a , k ⊥ a ) f ¯ q/p ↑ ( x b , k ⊥ b ) d ˆ σ q ¯ q → c ¯ c dM c ¯ c (cid:21)(cid:41) , (3)where q = u, d, s, ¯ u, ¯ d, ¯ s and the gluon and quark densities have been replaced by transversemomentum dependent gluon and quark PDFs. In Eq. 3, f g/p ↑ ( ↓ ) ( x, k ⊥ ) is the TMDPDFdescribing the distribution of gluons in proton which is transversely polarized w.r.t thebeam axis with the polarization being upwards (downwards) with respect to the productionplane. For a general value of the transverse spin S ⊥ , it is parametrised in terms of the gluonSivers function (GSF) ∆ N f g/p ↑ , as follows: f g/p ↑ ( x, k ⊥ , S ⊥ ; Q ) = f g/p ( x, k ⊥ ; Q ) − f ⊥ i T ( x, k ⊥ ; Q ) (cid:15) ab k a ⊥ S b ⊥ M p (4)= f g/p ( x, k ⊥ ; Q ) + 12 ∆ N f g/p ↑ ( x, k ⊥ , S ⊥ ; Q )= f g/p ( x, k ⊥ ; Q ) + 12 ∆ N f g/p ↑ ( x, k ⊥ ; Q ) (cid:15) ab k a ⊥ S b ⊥ k ⊥ . Any non-zero TSSA in the process considered would primarily arise due to an azimuthalanisotropy in the distribution of gluon transverse momenta in the polarized proton. Thisanisotropy is parametrised by the gluon Sivers distribution.Following Ref. [40] we can then write the numerator and denominator of Eq. 1 as,5 σ ↑ dyd q T − d σ ↓ dyd q T = F J/ψ s (cid:90) [ dM c ¯ c d k ⊥ a d k ⊥ b ] δ ( k ⊥ a + k ⊥ b − q T ) × (cid:40) ∆ N f g/p ↑ ( x a , k ⊥ a ) f g/p ( x b , k ⊥ b )ˆ σ gg → c ¯ c ( M c ¯ c )+ (cid:88) q (cid:20) ∆ N f q/p ↑ ( x a , k ⊥ a ) f ¯ q/p ( x b , k ⊥ b )ˆ σ q ¯ q → c ¯ c ( M c ¯ c ) (cid:21)(cid:41) (5)and d σ ↑ dyd q T + d σ ↓ dyd q T = 2 F J/ψ s (cid:90) [ dM c ¯ c d k ⊥ a d k ⊥ b ] δ ( k ⊥ a + k ⊥ b − q T ) × (cid:40) f g/p ↑ ( x a , k ⊥ a ) f g/p ( x b , k ⊥ b )ˆ σ gg → c ¯ c ( M c ¯ c )+ (cid:88) q (cid:20) f q/p ↑ ( x a , k ⊥ a ) f ¯ q/p ( x b , k ⊥ b )ˆ σ q ¯ q → c ¯ c ( M c ¯ c ) (cid:21)(cid:41) (6)with, x a,b = M c ¯ c √ s e ± y . (7)Here, y and q T are the rapidity and transverse momentum of the J/ψ and we consider theplane of production of the
J/ψ to be perpendicular to the proton spin S ⊥ . The partoniccross-sections for production of a c ¯ c pair of mass M are given by [41]ˆ σ gg → c ¯ c ( M ) = πα s s (cid:20)(cid:18) v + 116 v (cid:19) ln 1 + √ − v − √ − v − (cid:18)
74 + 3116 v (cid:19) √ − v (cid:21) (8)and ˆ σ q ¯ q → c ¯ c ( M ) = 29 (cid:18) πα s s (cid:19) (cid:18) v (cid:19) √ − v, (9)where v = m c M . B. Parametrization of the TMDs
For the predictions with the two DMP fits [33], we adopt the same functional forms forthe TMDs using which they were extracted. We use the standard form for the unpolarizedTMDPDF with a factorized Gaussian k ⊥ -dependence, f i/p ( x, k ⊥ ; Q ) = f i/p ( x, Q ) 1 π (cid:104) k ⊥ (cid:105) e − k ⊥ / (cid:104) k ⊥ (cid:105) (10)6ith (cid:104) k ⊥ (cid:105) = 0.25 GeV . The Sivers function is parameterized as [42]∆ N f i/p ↑ ( x, k ⊥ ; Q ) = 2 N i ( x ) f i/p ( x, Q ) h ( k ⊥ ) e − k ⊥ / (cid:104) k ⊥ (cid:105) π (cid:104) k ⊥ (cid:105) (11)with, N i ( x ) = N i x α i (1 − x ) β i ( α i + β i ) α i + β i α α i i β β i i (12)and h ( k ⊥ ) = √ e k ⊥ M e − k ⊥ /M , (13)where N i , α i , β i and M are all parameters determined by fits to data and e is Euler’snumber.As mentioned in the introduction, the two DMP extractions of the GSF, namely SIDIS1and SIDIS2 were obtained by fitting to data on TSSA in p ↑ p → π + X at RHIC with quarkSivers function (QSFs) extracted earlier from SIDIS data being used to account for thequark contribution to A N . In obtaining the SIDIS1 fits, only the u and d flavour were takeninto account, using the data on pion production from the HERMES experiment and positivehadron production from the COMPASS experiment. SIDIS2 parameters were obtained usingflavour segregated data on pion and kaon production so here all three light flavours weretaken into account. Furthermore, the QSFs used in the SIDIS1 fit were obtained with theset of fragmentation functions by Kretzer [43] and those of SIDIS2 were obtained with theset by de Florian, Sassot and Stratmann (DSS) [44]. The values of the parameters of thetwo fits are given in Table I. We give predictions for TSSA using these. SIDIS1 N g = 0 . α g = 2 . β g = 2 . ρ = 0 . (cid:104) k ⊥ (cid:105) = 0 . GeV SIDIS2 N g = 0 . α g = 0 . β g = 1 . ρ = 0 . ρ = M / ( (cid:104) k ⊥ (cid:105) + M ). C. TMD Evolution
The QCD evolution formalism of the unpolarised TMD and the Sivers function has beenobtained for DY and SIDIS [10] both of which involve a color singlet photon. One expectsthe TMD evolution of TMDs to be different for more complicated processes. However, since
J/ψ is a color singlet, in this case one can assume the TMD evolution of Ref. [10] for a7reliminary assessment of the effect of evolution on asymmetries. A more rigorous approachto TMD evolution for quarkonium will be closely related to the issue of validity of TMDfactorization for quarkonium which is not yet established.The energy evolution of a TMDPDF F ( x, k ⊥ ; Q ) is best described through its Fouriertransform into coordinate space which is given by F ( x, b ; Q ) = (cid:90) d k ⊥ e − i(cid:126)k ⊥ .(cid:126)b ⊥ F ( x, k ⊥ ; Q ) . (14)The evolution of b -space TMDPDFs can then be written as, F ( x, b ; Q f ) = F ( x, b, Q i ) R pert ( Q f , Q i , b ∗ ) R NP ( Q f , b ) , (15)where, R pert is the perturbatively calculable part of the evolution kernel, R NP is a nonpertur-bative Sudakov factor and b ∗ = b/ (cid:112) b/b max ) is used to stitch together the perturbativepart of the kernel, which is valid for b << b max , with the nonperturbative part, which isvalid for large b . Following Ref. [10], we choose an initial scale Q i = c/b ∗ for TMD evolution.Here c = 2 e − γ E where γ E is the Euler-Mascheroni constant.Setting Q i = c/b ∗ and Q f = Q , the perturbative part can be written as, R pert = exp (cid:26) − (cid:90) Qc/b ∗ dµµ (cid:18) A ln Q µ + B (cid:19)(cid:27) , (16)where A = Γ cusp and B = γ V are anomalous dimensions that can be expanded perturba-tively. The expansion coefficients with the appropriate gluon anomalous dimensions at NLLare [10] A (1) = C A , (17) A (2) = 12 C A (cid:18) C A (cid:18) − π (cid:19) − C A N f (cid:19) , (18) B (1) = − (cid:18) C A − N f (cid:19) . (19)The nonperturbative part of the evolution kernel, R NP is R NP = exp (cid:26) − b (cid:18) g TMD1 + g QQ (cid:19)(cid:27) , (20)where, g is a factor which takes the same value for all quark TMDPDFs [10, 11] and g isTMDPDF specific and is proportional to the intrinsic transverse momentum width of the8articular TMDPDF at the momentum scale Q . In case of gluon TMDPDFs, g is to bemultiplied by a factor of C A C F [45]. Assuming a factorized Gaussian form at scale Q , we have g unpol1 = (cid:104) k ⊥ (cid:105) Q g Sivers1 = (cid:104) k s ⊥ (cid:105) Q . (21)Expanding the TMDPDF F ( x, b ; Q ) at the initial scale in terms of its correspondingcollinear density at leading order (LO), for the unpolarized TMDPDF we get, f i/p ( x, b ; Q ) = f i/p ( x, c/b ∗ ) exp (cid:26) − (cid:90) Qc/b ∗ dµµ (cid:18) A ln Q µ + B (cid:19)(cid:27) (22) × exp (cid:26) − b (cid:18) g unpol1 + g QQ (cid:19)(cid:27) . In the case of the Sivers function, the evolution of its derivative in b -space can be writtenin the form of Eq. 15 leading to, f (cid:48)⊥ i T ( x, b ; Q ) = M p b T i,F ( x, x, c/b ∗ ) exp (cid:26) − (cid:90) Qc/b ∗ dµµ (cid:18) A ln Q µ + B (cid:19)(cid:27) (23) × exp (cid:26) − b (cid:18) g Sivers1 + g QQ (cid:19)(cid:27) , where, f (cid:48)⊥ T ( x, b ; µ ) ≡ ∂f ⊥ T ∂b . Here the Qiu-Sterman function for parton i , T i,F ( x, x, µ ) isobtained when expanding f (cid:48)⊥ T ( x, x, µ ) at LO and can be parametrized as T i,F ( x, x, µ ) = N i ( x ) f i/p ( x, µ ) (24)with N i ( x ) having the same form as in Eq. 12.The expressions for the TMDs in k ⊥ -space can be obtained by Fourier transforming the b -space expressions: f i/p ( x, k ⊥ ; Q ) = 12 π (cid:90) ∞ db bJ ( k ⊥ b ) f i/p ( x, b ; Q ) (25) f ⊥ i T ( x, k ⊥ ; Q ) = − πk ⊥ (cid:90) ∞ db bJ ( k ⊥ b ) f (cid:48)⊥ i T ( x, b ; Q ) . (26)The above expression for the Sivers function is related to ∆ N f i/p ↑ ( x, k ⊥ ; Q ) through Eq. 4.For the purpose of studying the effect of the transverse-momentum-dependent evolutionof the densities on the asymmetry predictions, we need to use the BV models for the GSFsince there are no available fits of the GSF that take into account TMD evolution. Wetherefore consider the following two models of the GSF wherein the x -dependent term N g ( x ),is modeled on that of the quarks: 9. BV (A): N g ( x ) = {N u ( x ) + N d ( x ) } /
22. BV (B): N g ( x ) = N d ( x )where both N u ( x ) and N d ( x ) are of the form given in Eq. 12. These models were first usedin Ref. [34] by Boer and Vogelsang and hence we will refer to these as the BV models.For the predictions with TMD evolved densities using the above two models, we use thefollowing set of parameters given in Ref. [10], N u = 0 . α u = 1 . β u = 4 . (cid:104) k ⊥ (cid:105) Q = 0 .
38 GeV N d = − . α d = 1 . β d = 4 . (cid:104) k s ⊥ (cid:105) Q = 0 .
282 GeV b max = 1 . − g = 0 .
16 GeV TABLE II: Parameters for TMD evolved densities from Ref. [10].
The predictions made with the TMD evolved densities using the BV models are comparedwith those obtained with DGLAP evolved densities using the same models with the followingset of parameters given in [46] N u = 0 . α u = 0 . β u = 2 . (cid:104) k ⊥ (cid:105) = 0 .
25 GeV N d = − . α d = 0 . β d = 0 . M = 0 . III. RESULTS
In this section, we present predictions of transverse single spin asymmetry in p + p ↑ → J/ψ + X , obtained using the DMP fits [33] and the BV models [34] of the GSF and cor-responding best fit parameters of QSFs. The QSFs corresponding to SIDIS-1 and SIDIS-2are given in ref. [20] and ref. [42] respectively. The best fits of QSF corresponding to BVmodels are given in ref. [46] for DGLAP evolved densities and in ref. [10] for the TMDevolved ones. Our predictions of TSSA are given for three different centre of mass energies √ s = 115 GeV (AFTER@LHC), 200 GeV (RHIC1) and 500 GeV (RHIC2). We presentasymmetry predictions as a function of: (i) the transverse momentum q T , with the rapid-ity integrated over − . ≤ y ≤ . √ s = 115 GeV, 2 ≤ y ≤ ≤ y ≤ . s = 200 GeV, and 2 ≤ y ≤ ≤ y ≤ √ s = 500 GeV (ii) the rapidity y , withthe transverse momentum integrated in the range 0 < q T < . ≤ η ≤
4. For convenience, we will referto (i) and (ii) as q T -asymmetry and y -asymmetry respectively. We then present a compar-ison of asymmetries estimated using DGLAP evolved densities with those obtained usingTMD-evolved densities in order to study the effect of TMD evolution. We then comparethe asymmetries obtained using the aforementioned fits and models, with measurements ofTSSA in p + p ↑ → J/ψ + X at √ s = 200 GeV, performed by the PHENIX collaboration atRHIC [36]. These measurements were performed in the forward (1 . ≤ y ≤ . − . ≤ y ≤ − .
2) and midrapidity ( − . ≤ y ≤ .
35) regions with 0 ≤ q T ≤ . . ≤ y ≤ .
0) that will become accessible under the proposed fsPHENIX upgrade [47, 48].
Gluon TotalSIDIS1SIDIS2BV (cid:72) A (cid:76) BV (cid:72) B (cid:76) (cid:45) (cid:45) q T (cid:72) GeV (cid:76) A N (cid:60) y (cid:60) (cid:61)
200 GeV (cid:45) (cid:45) q T (cid:72) GeV (cid:76) (cid:60) y (cid:60) (cid:61)
200 GeV
FIG. 1: Predictions of q T -asymmetry in p + p ↑ → J/ψ + X at RHIC1 ( √ s = 200 GeV) energy usingDGLAP evolved densities with BV(A), BV(B), DMP-SIDIS1 and DMP-SIDIS2 parameters. Leftpanel and right panel show q T -asymmetry integrated over the range 2 ≤ y ≤ ≤ y ≤ . To assess the contribution of the QSF over the GSF to the asymmetry, we have com-pared total asymmetry (contribution of both quarks and gluon Sivers functions) with thecontribution of gluon Sivers function to the asymmetry and found the contribution of QSF11
IDIS1SIDIS2BV (cid:72) A (cid:76) BV (cid:72) B (cid:76) (cid:45) (cid:45) (cid:45) (cid:45) y A N (cid:163) q T (cid:163) (cid:61)
200 GeV
FIG. 2: Predictions for y -asymmetry in p + p ↑ → J/ψ + X at RHIC1 ( √ s = 200 GeV) energyusing DGLAP evolved densities with BV(A), BV(B), DMP-SIDIS1 and DMP-SIDIS2 parameters. q T integration range is 0 ≤ q T ≤ . to be negligible in all cases. In Fig. 1, we show this comparison for DMP fits and the BVmodels of the GSF, at √ s = 200 GeV. We observe that the contribution of the QSF to theasymmetry is indeed very small as compared to contribution of GSF. This assures that theuse of this process as a probe of the gluon Sivers function will not be compromised. In theremaining figures, we have considered contribution from both QSF and GSF. In Fig. 2, weshow rapidity dependence of asymmetry predictions obtained using different sets of DGLAPevolved densities, i.e., the DMP fits and the BV models of the GSF, at √ s = 200 GeV. Weobserve that the signs of asymmetry obtained using BV parameters and more recent directlyfitted DMP parameters are opposite. This is expected as in BV models, the gluon Siversfunction is modelled after quark Sivers function and the d -quark Sivers parameters have anegative sign as shown in Tables II and III. Further, the magnitude of asymmetry obtainedusing the DMP fits is smaller than that obtained using the BV models. Of the two DMPfits, SIDIS1 gives the larger asymmetry estimates with peak values of about 35% for the q T -asymmetry (with the rapidity range 3 ≤ y ≤ . y -asymmetry.SIDIS2 on the other hand gives much smaller asymmetries with peak values of about 2%for the q T -asymmetry (with the rapidity range 3 ≤ y ≤ .
8) and 2% for the y -asymmetry.This large difference in the peak magnitudes of the asymmetry between SIDIS1 and SIDIS2fits can be understood by looking at the x a -region which contributes to the peaks. Forboth SIDIS1 and SIDIS2, the peak occurs for y > x region12 a > ∼ . y -asymmetry scales with P L / √ s . It should be noted that in Fig. 3b and4b, the y -asymmetry peaks in negative y region for AFTER@LHC c.m energy. This isdue to the fact that AFTER@LHC is a fixed target experiment and we have taken y cms to be positive in the (unpolarized) beam direction. This is in contrast to RHIC1 andRHIC2 curves, where we have used the convention followed by PHENIX experiment whererapidity is considered to be positive in the forward hemisphere of the polarized proton.The scaling of y -asymmetry with P L / √ s is because the collinear PDFs mostly cancel outbetween the numerator and the denominator and the y -dependence of the asymmetry ismostly determined by the remaining factor N g ( x ), with x a having a direct correspondencewith y (c.f. Eq. 7). This cancellation is of course not absolute, as the integration over theinvariant mass M c ¯ c dilutes the correspondence of x a with the rapidity and allows the collinearPDFs to affect the y -dependence, but we have verified that this effect is small. It must alsobe mentioned that the assumption that the k ⊥ dependence of the TMD is factorized, helpsthis cancellation as the integrals over the transverse momenta of the partons do not dependon √ s and simply produce an overall constant independent of both y and √ s . We find peakasymmetry values of about 27% with SIDIS1 and 1 .
9% with SIDIS2.In the case of the q T -asymmetries, shown in Fig. 3a and 4a, we find that the functionalform of the q T dependence remains the same up to an overall factor that depends on √ s and the rapidity range. This is also a reflection of the factorized k T -dependence that wehave assumed for the TMDPDFs. With SIDIS1, we find that peak asymmetry occurs at q T ∼ . . √ s = 115 (with − . ≤ y ≤ .
2) , 200 (with 3 ≤ y ≤ .
8) and 500 (with 3 ≤ y ≤
4) GeV respectively, whilefor SIDIS2 we get substantially lower asymmetries with peak values of 1.1%, 2.2% and 1.9%for √ s = 115 (with − . ≤ y ≤ . ≤ y ≤ .
8) and 500 (with 3 ≤ y ≤
4) GeVrespectively at q T ∼ . . (cid:61)
115 GeV, (cid:45) (cid:60) y (cid:60) s (cid:61)
200 GeV, 2 (cid:60) y (cid:60) (cid:61)
200 GeV, 3 (cid:60) y (cid:60) (cid:61)
500 GeV, 2 (cid:60) y (cid:60) (cid:61)
500 GeV, 3 (cid:60) y (cid:60) q T (cid:72) GeV (cid:76) A N SIDIS1 s (cid:61)
115 GeVs (cid:61)
200 GeVs (cid:61)
500 GeV (cid:45) (cid:45) y A N SIDIS10 (cid:60) q T (cid:60) FIG. 3: Predictions for asymmetry as a function of q T (left panel) and y (right panel) obtained usingthe DMP-SIDIS1 GSF [33] parameters for all the three centre of mass values considered ( √ s = 115GeV, 200 GeV, 500 GeV). Integration ranges for rapidity in left panel are − . ≤ y ≤ . √ s = 115 GeV, 2 ≤ y ≤ ≤ y ≤ . √ s = 200 GeV and 2 ≤ y ≤ ≤ y ≤ √ s = 500 GeV. In right panel, asymmetry predictions obtained using the DMP-SIDIS1 GSF [33]as a function of y are presented for all three centre of mass values ( √ s = 115 GeV, 200 GeV, 500GeV). Integrati on range is for q T is 0 ≤ q T ≤ . y region for AFTER@LHC energy as we have used the convention for fixed target experiments asexplained in the Section III) evolved densities. We do this using the BV (A) and BV (B) models of the GSF (c.f. SectionII C). We find that the inclusion of TMD evolution causes the asymmetry predictions tosubstantially decrease for both models. Furthermore the P L / √ s scaling of the y -asymmetryis also affected by TMD evolution as can be seen in Fig. 6, where we show the asymmetriesobtained with the TMD evolved BV (B) model, for all three c.o.m energies. While thepeak of the asymmetry does shift to larger rapidities with increasing √ s , the magnitude ofthe peak varies. This is due to the fact that the k T -dependence of the TMDs is no morefactorized, but is instead affected by the x -dependence of the collinear PDFs through the c/b ∗ prescription (c.f., Eq. 23, 24).In Figs. 7 and 8, we compare our predictions with the asymmetry measured at thePHENIX experiment [36]. In Fig. 7, we compare the asymmetry predictions obtained usingthe DMP fits with the data and find that they lie well within the uncertainties. In the forwardregion, SIDIS1 and SIDIS2 give an asymmetries of about 1 .
2% and 0 .
9% respectively. In14
15 GeV, (cid:45) (cid:60) y (cid:60) (cid:60) y (cid:60) (cid:60) y (cid:60) (cid:60) y (cid:60) (cid:60) y (cid:60) q T (cid:72) GeV (cid:76) A N SIDIS2
115 GeV200 GeV500 GeV (cid:45) (cid:45) y A N SIDIS20 (cid:60) q T (cid:60) FIG. 4: Predictions for asymmetry as a function of q T (left panel) and y (right panel) obtained usingthe DMP-SIDIS2 GSF [33] parameters for all three centre of mass values considered ( √ s = 115GeV, 200 GeV, 500 GeV). Integration ranges for rapidity in left panel are − . ≤ y ≤ . √ s = 115 GeV, 2 ≤ y ≤ ≤ y ≤ . √ s = 200 GeV and 2 ≤ y ≤ ≤ y ≤ √ s = 500 GeV. In right panel, asymmetry predictions obtained using the DMP-SIDIS2 GSF [33]are plotted as a function of y for all three centre of mass values considered ( √ s = 115 GeV, 200GeV, 500 GeV). Integration range for q T is 0 ≤ q T ≤ . y region for AFTER@LHC energy as we have used the convention for fixed targetexperiments as explained in the Section III). Fig. 8, we do the same using the DGLAP and TMD evolved BV models. With the DGLAPevolved models, BV (A) gives asymmetries which lie within the uncertainties, whereas BV(B) gives an asymmetry well outside the uncertainties for the forward region. However, theasymmetry predictions obtained with the TMD evolved BV models are substantially smallerthan those given by DGLAP evolved BV models and are negligible in all rapidity regions.Thus, although the PHENIX results for the asymmetry are compatible with zero, theerrors are still large and allow for percentage level asymmetries as given by the DMP fits.Furthermore, they cover a very limited kinematic range and do not rule out larger asym-metries in more forward regions. The proposed fsPHENIX upgrade [47, 48] will expandthe forward coverage of the detecter to the region η ≤
4. With this in mind, in Fig. 9,we show asymmetry predictions obtained from the DMP fits for the expanded rapidity re-gion that will be covered by the upgrade. The expected statistical error for each point,which is given by δA N = (cid:112) (1 − A N ) /N , was calculated assuming 1 pb − of data. Here15 V (cid:72) A (cid:76) (cid:45) DGLAPBV (cid:72) B (cid:76) (cid:45) DGLAPBV (cid:72) A (cid:76) (cid:45) TMDBV (cid:72) B (cid:76) (cid:45) TMD (cid:45) (cid:45) q T (cid:72) GeV (cid:76) A N (cid:60) y (cid:60) (cid:61)
200 GeV BV (cid:72) A (cid:76) (cid:45) DGLAPBV (cid:72) B (cid:76) (cid:45) DGLAPBV (cid:72) A (cid:76) (cid:45) TMDBV (cid:72) B (cid:76) (cid:45) TMD (cid:45) (cid:45) (cid:45) (cid:45) y (cid:60) q T (cid:60) (cid:61)
200 GeV
FIG. 5: Comparison of asymmetry predictions obtained using TMD evolved BV models [34] withthose obtained using DGLAP evolved BV models. Left panel shows the q T -asymmetry integratedover 2 ≤ y ≤ y -asymmetry integrated over 0 ≤ q T ≤ . s (cid:61)
115 GeVs (cid:61)
200 GeVs (cid:61)
500 GeV (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) y A N BV (cid:72) B (cid:76) (cid:45) TMD0 (cid:60) q T (cid:60) FIG. 6: Predictions for y -asymmetry obtained using the TMD evolved BV (B) model of GSF [34] forall three centre of mass values considered ( √ s = 115 GeV, 200 GeV, 500 GeV) with q T integrationrange 0 ≤ q T ≤ . N = σ pp → J/ψ × B.R(
J/ψ → µ + µ − ) × F µ + µ − J/ψ × L , indicates the number
J/ψ ’s that are de-tected. L is the integrated luminosity, which we choose to be 1 pb − here and F µ + µ − J/ψ is thegeometric factor accounting for the planned detector acceptance of leptons: − . ≤ η l ≤ . σ pp → J/ψ was calculated using the colour evaporation model (CEM), nor-malised to the total cross-section given in [49].In Fig. 9, we would like to highlight the widely differing behavior of the SIDIS1 and SIDIS2asymmetries with respect to the choice of the rapidity cuts. Note that use of these fits for16 (cid:243) (cid:243)(cid:231) (cid:231)(cid:237) (cid:237) (cid:230) (cid:230) (cid:230)(cid:224) (cid:224) (cid:224) (cid:237)
PHENIX 2006 (cid:43) (cid:231)
PHENIX 2008 (cid:243)
PHENIX 2006 (cid:230)
SIDIS (cid:45) (cid:224) SIDIS (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) x F A N FIG. 7: Comparison of PHENIX measurements [36] of TSSA in p + p ↑ → J/ψ + X with predictionsobtained using the DMP fits, SIDIS1 and SIDIS2 [33]. The points for the combined (2006 +2008)data have been offset by 0.01 in x F for visibility. Asymmetry measurements are in the forward(1.2 < y < < y < -1.2) and midrapidity ( | y | < .
35) regions with 0 ≤ q T ≤ . GSF to make predictions for our process assumes TMD factorisation and universality ofGSF. These have yet to be established . Even then these two can be used for demonstrationpurposes, as examples of two possible GSF, which have widely different x -dependencies.Taking these as examples we demonstrate how a study of asymmetry with different choicesof rapidity cuts in the forward region, can probe the x -dependence of the GSF. When thewhole accessible region − . ≤ y ≤ . σ of the statistical error. For the region 1 . ≤ y ≤ . ≤ y ≤
3, SIDIS1 gives anasymmetry of about 7%, which is almost five times of the prediction given by SIDIS2. Forthe region 3 ≤ y ≤ .
8, the difference is even larger with SIDIS1 giving as asymmetry of 22%,which greater than that of SIDIS2 by a factor of 13. This difference is due to the different x -dependence of the fits. The region y > x ↑ > ∼ .
11 where SIDIS1 ismuch larger than SIDIS2. This difference however, is not seen in the asymmetry withoutrapidity cuts, i.e., considering the whole region − . ≤ y ≤ .
8, since the cross-section17 (cid:243) (cid:243)(cid:231) (cid:231)(cid:237) (cid:237)(cid:230) (cid:230) (cid:230)(cid:224) (cid:224) (cid:224)(cid:237)
PHENIX 2006 (cid:43) (cid:231)
PHENIX 2008 (cid:243)
PHENIX 2006 (cid:230) BV (cid:72) A (cid:76) (cid:224) BV (cid:72) B (cid:76) (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) x F A N (cid:243) (cid:243) (cid:243)(cid:231) (cid:231)(cid:237) (cid:237)(cid:230) (cid:230) (cid:230)(cid:224) (cid:224) (cid:224)(cid:237) PHENIX 2006 (cid:43) (cid:231)
PHENIX 2008 (cid:243)
PHENIX 2006 (cid:230) BV (cid:72) A (cid:76) (cid:45) TMD (cid:224) BV (cid:72) B (cid:76) (cid:45) TMD (cid:45) (cid:45) (cid:45) (cid:45) (cid:45) x F FIG. 8: Comparison of PHENIX measurements [36] of TSSA in p + p ↑ → J/ψ + X with predictionsobtained using BV models of the GSF [34]. The points for the combined (2006 +2008) data havebeen offset by 0.01 in x F for visibility. Left panel shows comparison with DGLAP evolved BVparameters and right panel shows the same for TMD evolved BV parameters. drops rapidly in the large rapidity region. We therefore see that a study of the dependenceof the measured asymmetries on the rapidity region over which measurement is made, cangive insight into the x -dependence of the gluon Sivers function. The numerical values of thecross-section, the asymmetries and their associated error are given in Table IV. (cid:224)(cid:224) (cid:230)(cid:230)(cid:236)(cid:236) (cid:225)(cid:225) (cid:231)(cid:231) (cid:243)(cid:243)(cid:237)(cid:237) (cid:37) (cid:37) (cid:37) (cid:37) x F A N (cid:242)(cid:242) (cid:37) (cid:37) (cid:37) (cid:224) : 1.2 (cid:163) y (cid:163) (cid:230) : 2.0 (cid:163) y (cid:163) (cid:242) : 3.0 (cid:163) y (cid:163) (cid:236) : (cid:45) (cid:163) y (cid:163) (cid:236) : SIDIS1 (cid:237) : SIDIS2 s (cid:61)
200 GeV0 (cid:163) q T (cid:163) FIG. 9: Predictions for asymmetry in forward region that would be accessible with the fsPHENIXupgrade [47, 48]. Error bars indicate expected statistical errors calculated assuming 1 pb − of data. egion (cid:104) x F (cid:105) x ↑ region σ J/ψ × B.R( µ + µ − ) (nb) δA N A SIDIS1 N A SIDIS2 N − . < y < . . < y < . < y < < y < . √ s = 200 GeV), along withexpected statistical error (assuming 1 pb − of data) and approximate x -region probed. All numbersgiven with cuts on lepton rapidity: − . < y l < .
0. Transverse momentum region: 0 ≤ q T ≤ . Another experiment that might help study the GSF in pp ↑ → J/ψ + X over kinematicregions not covered by the PHENIX study would be the proposed fixed target experimentAFTER@LHC. This would have a high enough luminosity to make precise measurementsof the asymmetry [50–55]. Such a fixed target experiment would have a centre of massenergy √ s = 115 GeV with an integrated luminosity of up to 20 fb − with one year ofdata taking. In such an experiment the pp centre of mass would be moving with respectto the lab frame with a rapidity y cms = 4 .
8, allowing large x regions of the target to beprobed with the coverage of the ALICE or LHCb detectors. A polarized target would,therefore, offer the possibility of probing the large x ↑ ( > ∼ .
3) region where the DMP fits donot constrain the GSF.
J/ψ production rates obtained with the leading order (LO) CEMindicate that, with an integrated luminosity of 20 fb − , it should be possible to measure theasymmetry with permille precision in the low- q T region with rapidity range − . ≤ y ≤ − . . ≤ y cms ≤ . . < ∼ x ↑ < ∼ . IV. CONCLUSIONS
In this paper, we have presented predictions of TSSA in p + p ↑ → J/ψ + X at theRHIC centre of mass energy √ s = 200 GeV, at which TSSA in J/ψ production has beenmeasured by the PHENIX experiment as well as at √ s = 115 GeV, corresponding to theproposed polarized scattering experiment AFTER@LHC and √ s = 500 GeV, correspondingto the higher RHIC energy using two different parameterizations of gluon Sivers functionand in different experimentally accessible kinematic regions. Measurement of TSSA in J/ψ p ↑ p → π + X , referred to here as the DMP parame-ter sets [33]. These extractions of GSF were obtained without taking into account TMDevolution. Hence, we have not used these parameters of GSF to assess the effect of TMDevolution. For the comparison of asymmetries calculated using DGLAP evolved and TMDevolved densities, we have used earlier models of the GSF (referred to as BV models in thiswork), which express GSF in terms of u and d quark Sivers function and were used in ourprevious work.Our results show that the asymmetry obtained using BV parameterization is negativewhereas asymmetry obtained using SIDIS parameterization is positive. This can be un-derstood considering the fact that d quark Sivers function has negative sign and in BVmodel, GSF is written in terms of QSF. The q T distribution of asymmetry we obtainedwith BV parameterization is large in magnitude as compared to asymmetry obtained us-ing DMP parameterizations. As far as the y distribution of asymmetry is concerned, thepeak magnitudes of asymmetry are similar for all c.o.m energies, with the peak shiftingtowards larger rapidity values with increasing c.o.m energy. Comparison of our predictionsof asymmetry with PHENIX measurement gives interesting results. When DGLAP evolvedTMDs are used, our results with DMP-SIDIS1, DMP-SIDIS2 and BV (A) parameterizationsare in agreement with the asymmetry data in forward, backward and midrapidity regionswhereas the results obtained using BV (B) parameterization are not within experimentaluncertainties of data points. However, when effect of TMD evolution is taken into account,results obtained using BV (B) parameterization also fall within experimental uncertainties.It may be worthwhile to obtain fits of GSF taking into account TMD evolution and comparepredictions of asymmetry obtained using those with the predictions presented here.As mentioned earlier, since the uncertainties in data points are large, more data will beneeded to constrain the GSF. The predictions made in this work are based on the first everdirectly fitted parameters of gluon Sivers function and assume a generalized factorizationexpression within colour evaporation model of charmonium production. A more detailedanalysis investigating the dependence of our results on charmonium production mechanism,20hich is still an open question, is under study and will be reported in future. Apart fromuncertainties arising from underlying assumptions which include TMD factorization forquarkonium production, which has not yet been established, universality of GSF and choiceof a particular production mechanism, another issue that is unresolved so far, there may befurther limitations due to restricted region of validity of the parameter sets used. However,these studies to understand the GSF and resulting TSSA in pp ↑ along with studies probingGSF in open flavour production are expected to play a crucial role in constraining the gluonspin dynamics. V. ACKNOWLEDGEMENTS
We would like to thank J.P. Lansberg for his very useful comments and suggestions.A.M. and B.S. would like to thank DST, India for financial support under the projectno.EMR/2014/0000486 and UGC-BSR under F.7-130/2007(BSR). A.M. would also like tothank the Theoretical Physics Department, CERN, Geneva, where part of this work wasdone, for their kind hospitality. R.M.G. wishes to acknowledge support from the Departmentof Science and Technology, India under Grant No. SR/S2/JCB-64/2007 under the J.C.Bose Fellowship scheme. A.K. would like to thank the Department of Physics, University ofMumbai, for their kind hospitality. [1] A. Adare et al. (PHENIX), Phys. Rev.
D90 , 012006 (2014), 1312.1995.[2] R. D. Klem, J. E. Bowers, H. W. Courant, H. Kagan, M. L. Marshak, E. A. Peterson, K. Rud-dick, W. H. Dragoset, and J. B. Roberts, Phys. Rev. Lett. , 929 (1976).[3] J. Antille, L. Dick, L. Madansky, D. Perret-Gallix, M. Werlen, A. Gonidec, K. Kuroda, andP. Kyberd, Phys. Lett. B94 , 523 (1980).[4] G. L. Kane, J. Pumplin, and W. Repko, Phys. Rev. Lett. , 1689 (1978).[5] A. Airapetian et al. (HERMES), Phys. Rev. Lett. , 152002 (2009), 0906.3918.[6] M. Alekseev et al. (COMPASS), Phys. Lett. B673 , 127 (2009), 0802.2160.[7] C. Adolph et al. (COMPASS), Phys. Lett.
B717 , 383 (2012), 1205.5122.
8] X. Qian et al. (Jefferson Lab Hall A), Phys. Rev. Lett. , 072003 (2011), 1106.0363.[9] U. D’Alesio and F. Murgia, Phys. Rev.
D70 , 074009 (2004), hep-ph/0408092.[10] M. G. Echevarria, A. Idilbi, Z.-B. Kang, and I. Vitev, Phys. Rev.
D89 , 074013 (2014),1401.5078.[11] J. Collins,
Foundations of perturbative QCD (Cambridge University Press, 2013),ISBN 9781107645257, 9781107645257, 9780521855334, 9781139097826, URL .[12] X.-d. Ji, J.-P. Ma, and F. Yuan, Phys. Lett.
B597 , 299 (2004), hep-ph/0405085.[13] X.-d. Ji, J.-p. Ma, and F. Yuan, Phys. Rev.
D71 , 034005 (2005), hep-ph/0404183.[14] A. Bacchetta, D. Boer, M. Diehl, and P. J. Mulders, JHEP , 023 (2008), 0803.0227.[15] D. Boer et al. (2011), 1108.1713.[16] A. Bacchetta and M. Radici, Phys. Rev. Lett. , 212001 (2011), 1107.5755.[17] D. W. Sivers, Phys. Rev. D41 , 83 (1990).[18] D. W. Sivers, Phys. Rev.
D43 , 261 (1991).[19] M. Anselmino, M. Boglione, U. D’Alesio, S. Melis, F. Murgia, and A. Prokudin, J. Phys. Conf.Ser. , 012062 (2011), 1012.3565.[20] M. Anselmino, M. Boglione, U. D’Alesio, A. Kotzinian, F. Murgia, and A. Prokudin, Phys.Rev.
D72 , 094007 (2005), [Erratum: Phys. Rev.D72,099903(2005)], hep-ph/0507181.[21] A. Bacchetta, F. Conti, and M. Radici, Phys. Rev.
D78 , 074010 (2008), 0807.0323.[22] L. P. Gamberg, G. R. Goldstein, and M. Schlegel, Phys. Rev.
D77 , 094016 (2008), 0708.0324.[23] I. Schmidt, J. Soffer, and J.-J. Yang, Phys. Lett.
B612 , 258 (2005), hep-ph/0503127.[24] M. Anselmino, M. Boglione, U. D’Alesio, E. Leader, and F. Murgia, Phys. Rev.
D70 , 074025(2004), hep-ph/0407100.[25] Z.-B. Kang and J.-W. Qiu, Phys. Rev.
D78 , 034005 (2008), 0806.1970.[26] A. Bacchetta, C. Bomhof, U. D’Alesio, P. J. Mulders, and F. Murgia, Phys. Rev. Lett. ,212002 (2007), hep-ph/0703153.[27] A. Schafer and J. Zhou, Phys. Rev. D88 , 014008 (2013), 1302.4600.[28] F. Yuan, Phys. Rev.
D78 , 014024 (2008), 0801.4357.[29] S. J. Brodsky, D. S. Hwang, and I. Schmidt, Phys. Lett.
B530 , 99 (2002), hep-ph/0201296.[30] S. J. Brodsky, D. S. Hwang, and I. Schmidt, Nucl. Phys.
B642 , 344 (2002), hep-ph/0206259.[31] R. M. Godbole, A. Kaushik, A. Misra, V. Rawoot, and B. Sonawane, Few Body Syst. , 96 D96 , 036011 (2017), 1705.04169.[33] U. D’Alesio, F. Murgia, and C. Pisano, JHEP , 119 (2015), 1506.03078.[34] D. Boer and W. Vogelsang, Phys. Rev. D69 , 094025 (2004), hep-ph/0312320.[35] R. M. Godbole, A. Kaushik, and A. Misra, Phys. Rev.
D94 , 114022 (2016), 1606.01818.[36] A. Adare et al. (PHENIX), Phys. Rev.
D82 , 112008 (2010), [Erratum: Phys.Rev.D86,099904(2012)], 1009.4864.[37] M. B. Gay Ducati and C. Brenner Mariotto, Phys. Lett.
B464 , 286 (1999), hep-ph/9908407.[38] F. Halzen and S. Matsuda, Phys. Rev.
D17 , 1344 (1978).[39] U. D’Alesio and F. Murgia, Prog. Part. Nucl. Phys. , 394 (2008), 0712.4328.[40] R. M. Godbole, A. Misra, A. Mukherjee, and V. S. Rawoot, Phys. Rev. D85 , 094013 (2012),1201.1066.[41] M. Gluck and E. Reya, Phys. Lett. , 453 (1978).[42] M. Anselmino, M. Boglione, U. D’Alesio, A. Kotzinian, S. Melis, F. Murgia, A. Prokudin, andC. Turk, Eur. Phys. J.
A39 , 89 (2009), 0805.2677.[43] S. Kretzer, Phys. Rev.
D62 , 054001 (2000), hep-ph/0003177.[44] D. de Florian, R. Sassot, and M. Stratmann, Phys. Rev.
D75 , 114010 (2007), hep-ph/0703242.[45] Z.-B. Kang, X. Liu, F. Ringer, and H. Xing (2017), 1705.08443.[46] M. Anselmino, M. Boglione, U. D’Alesio, S. Melis, F. Murgia, and A. Prokudin, in (2011), 1107.4446, URL .[47] K. N. Barish (PHENIX), J. Phys. Conf. Ser. , 012033 (2012).[48] E.-C. Aschenauer et al. (2015), 1501.01220.[49] A. Adare et al. (PHENIX), Phys. Rev. Lett. , 232002 (2007), hep-ex/0611020.[50] J.-P. Lansberg et al., PoS PSTP2015 , 042 (2016), 1602.06857.[51] M. Anselmino, U. D’Alesio, and S. Melis, Adv. High Energy Phys. , 475040 (2015),1504.03791.[52] S. J. Brodsky, F. Fleuret, C. Hadjidakis, and J. P. Lansberg, Phys. Rept. , 239 (2013),1202.6585.[53] D. Kikoa, M. G. Echevarria, C. Hadjidakis, J.-P. Lansberg, C. Lorc, L. Massacrier, C. M. uintans, A. Signori, and B. Trzeciak (2017), 1702.01546.[54] J. P. Lansberg, S. J. Brodsky, F. Fleuret, and C. Hadjidakis, Few Body Syst. , 11 (2012),1204.5793.[55] A. Rakotozafindrabe et al., Phys. Part. Nucl. , 336 (2014), 1301.5739., 336 (2014), 1301.5739.