Production of ρ^{0} meson with large p_T at NLO in heavy-ion collisions
aa r X i v : . [ nu c l - t h ] A ug Production of ρ meson with large p T at NLO in heavy-ion collisions Wei Dai, Ben-Wei Zhang ∗ , and Enke Wang School of Mathematics and Physics, China University of Geosciences (Wuhan), Wuhan 430074, China Key Laboratory of Quark & Lepton Physics (MOE) and Institute of Particle Physics,Central China Normal University, Wuhan 430079, China (Dated: October 14, 2018)Production of large transverse momentum ρ meson in high-energy nuclear collisions is inves-tigated for the first time at the next-leading-order in the QCD improved parton model. The ρ fragmentation functions (FFs) in vacuum at any scale Q are obtained, by evolving a newly devel-oped initial parametrization of ρ FFs at a scale Q = 1 . from a broken SU(3) model throughNLO DGLAP equations. The numerical simulations of p T spectra of ρ meson in the elementaryp + p collisions at NLO give a decent description of STAR p + p data. In A + A reactions the jetquenching effect is taken into account with the higher-twist approach by the medium-modified par-ton FFs due to gluon radiation in the quark-gluon plasma, whose space-time evolution is describedby a (3+1D) hydrodynamical model. The nuclear modification factors for ρ meson and its doubleratio with π ± nuclear modification in central Au + Au collisions at the RHIC are calculated andfound to be in good agreement with STAR measurement. Predictions of ρ nuclear modificationand the yield ratio ρ /π in central Pb+Pb at the LHC are also presented. It is shown that theratio ρ /π in central Pb+Pb will approach to that in p+p reactions when p T >
12 GeV.
PACS numbers: 12.38.Mh; 25.75.-q; 13.85.Ni
A new state of matter of deconfined quarks and gluons,the so-called quark-gluon plasma (QGP), is expected tobe created in heavy ion collisions (HIC) at very high col-liding energies. To study the creation and properties ofthe QGP, the jet quenching has been proposed, whichstates that when an energetic parton traveling throughthe hot/dense QCD medium, a substantial fraction of itsenergy should be losted and could in turn be used toobtain the temperature and density information of theQGP [1, 2]. Even though rapid developments of exper-iments and theories on new jet quenching observables,such as di-hadron [3, 4], photon triggered hadron [5, 6]and full jet observable [7–12], have emerged in the lastdecade, the suppression of inclusive hadron production,as the most intensively studied observable on jet quench-ing, is still indispensable to unravel the properties ofthe QCD medium. Recently, by comparing the theoreti-cal calculation with the measurements of the productionspectra and its suppression of π mesons which are themost commonly observed hadrons, the jet transport co-efficient ˆ q has been extracted to characterize the localproperties of the QCD medium probed by the energeticparton jets [13]. The higher twist multiple scattering ofthe jet quenching incorporated with perturbative quan-tum chromodynamics (pQCD) improved parton modelhas been developed and successfully described the π and η productions and their suppressions in A + A colli-sions [14–17].The study of the identified hadron spectra at high p T other than π and η in HIC can further constrain and castinsight into the hadron suppression pattern. Whereas a ∗ [email protected] relatively large amount of data on the yields of identifiedhadrons at large p T has been accumulated at the RHICand the LHC [18–20], there are still very few theoreticalstudies of hadrons with different types. An interestingtype of identified hadrons with available data is ρ me-son, which is heavier than π and η , and also consists ofthe similar constituent quarks. We notice that even thetheoretical calculations of the ρ productions in p+p col-lisions with large p T at both the RHIC and the LHC areabsent due to the lack of knowledge of parton fragmen-tation functions (FFs) for ρ in vacuum. In a previousstudy [16] we have paved the way to understand identi-fied hadron suppression pattern by calculating the pro-ductions of η meson and investigating the hadron yieldratios [16]. In this manuscript, we extend this study to ρ meson productions and the yield ratios of ρ and π inA+A collisions at the RHIC and the LHC. It is of greatinterest to see how the alteration of the jet chemistrybrought by the jet quenching will eventually affect the ρ production spectrum and the ratio of hadron yields [21–23].In this paper, firstly we employ a newly developed ini-tial parametrization of ρ FFs in vacuum at a startingscale Q = 1 . , which is provided by the SU (3)model of FFs of vector mesons [24, 25]. By evolving themthrough DGLAP evolution equations at NLO [26], we ob-tain parton FFs of ρ meson at any hard scale Q . Thetheoretical results of ρ productions in p+p collisions areprovided up to the next-to-leading order(NLO) in pQCDimproved parton model, and we find that they describethe experimental data rather well. Then we study ρ production in A + A collisions at both RHIC and LHCby including parton energy loss in the hot/dense QCDmedium in the framework of higher twist approach ofjet quenching [27–29]. In this approach, the energy lossdue to the multiple scattering suffered by an energeticparton traversing the medium are taken into account bytwist-4 processes, and the vacuum fragmentation func-tions are modified effectively in high-energy nuclear col-lisions. Therefore, we can compute numerically for thefirst time ρ meson yields in A + A collisions. We givea description of ρ nuclear modification factor R AA ( ρ )at large p T in Au + Au collisions at the RHIC to con-front against the experimental data by STAR Collabo-ration, and R AA ( ρ ) in Pb + Pb collisions at the LHCto give a theoretical prediction. The double ratio of R AA ( ρ ) /R AA ( π ± ) is calculated and found to be in goodagreement with the experimental data. Lastly we explorethe features of the ρ /π ratios in both p+p and A+Acollisions.In NLO pQCD calculation, the single hadron produc-tion can be factorized as the convolution of elementarypartonic scattering cross sections up to α , parton dis-tribution functions (PDFs) inside the incoming particlesand parton FFs to the final state hadrons [30]. We canexpress the formula symbolically as:1 p T dσ h dp T = Z F q ( p T z h ) · D q → h ( z h , p T ) dz h z h + Z F g ( p T z h ) · D g → h ( z h , p T ) dz h z h . (1)The above equation implies that the hadron yield in p+pcollision will be determined by two factors: the initial(parton-)jet spectrum F q,g ( p T ) and the parton fragmen-tation functions D q,g → h ( z h , p T ). In the following calcu-lations, we utilize CTEQ6M parametrization for protonPDFs [31], which has been convoluted with the elemen-tary partonic scattering cross sections up to α to obtain F q,g ( p T z h ). Here D q,g → h ( z h , p T ) represents the vacuumparton FFs, which denote the possibilities of scatteredquark or gluon fragmenting into hadron h with momen-tum fraction z h . They can be given by correspondingparametrization for different final-state hadrons. So po-tentially, we could predict all the identified hadron pro-ductions in p + p collision as long as the fragmentationfunctions are available. Note that the factorization scale,renormalization scale and fragmentation scale are usu-ally chosen to be the same and proportional to p T of theleading hadron in the final-state.To accurately determine the p + p reference, partonFFs in vacuum as a non-perturbative input, should beavailable. So far it is still impossible to derive parton FFsfrom the first-principle of QCD and a common practice isto make phenomenological parametrizations by compar-ing perturbative QCD calculations with the data. Un-like π and charged hadrons, until now there are veryfew satisfatory parametrizations of parton FFs for thevector mesons due to the paucity of the relevant data.Fortunately, a broken SU (3) model is recently proposedto provide a systematic description of the vector mesonsproduction [24, 25]. To reduced the complexity of the meson octet fragmentation functions, the SU (3) flavorsymmetry is introduced with a symmetry breaking pa-rameter. In addition, isospin and charge conjugation in-variance of the vector mesons ρ ( ρ + , ρ − , ρ ) are assumedto further reduce independent unknown quark FFs intofunctions named valence(V) and sea( γ ). The inputs ofvalence V ( x, Q ), sea γ ( x, Q ) and gluon D g ( x, Q ) FFsare parameterized into a standard polynomial at a start-ing low energy scale of Q = 1 . such as: F i ( x ) = a i x b i (1 − x ) c i (1 + d i x + e i x ) (2)These parameters are systematically fixed by fitting thecross section at NLO with the measurements of LEP( ρ , ω )and SLD( φ , K ⋆ ) at √ s = 91 . ρ FFs in vacuum at Q = 1 . arelisted and we obtain ρ FFs at any hard scale D q,g ( x, Q ) Q > ρ FFs D q,g ( x, Q ) are used in ournumerical simulations.We have plotted the parton FFs as functions of frag-menting fraction z h in the left panel of Fig. 1 at fixedscale of Q = 100 GeV , and also the parton FFs asfunctions of final state p T at fixed fragmenting fraction z h = 0 . ρ FFs decrease with z h , and FF of up quarkis much larger than that of strange quark, especially atlarge z region. At a typical value with z h = 0 .
6, we noticethat ρ FFs show a rather weak p T dependence. h z q , g D -7 -6 -5 -4 -3 -2 -1 u quarks quarkgluon =100 GeV Q (GeV) T p − − q , g D u quarks quarkgluon =0.6 h z FIG. 1: Left: parton FFs as functions of z h at fixed scale Q = 100 GeV ; Right: parton FFs as functions of p T at fixed z h = 0 . The existence of the ρ meson FFs at NLO allows usto calculate the inclusive vector meson productions as afunction of the final state hadron p T in pQCD at theaccuracy of NLO. Fig. 2 shows the confrontation of thetheoretical calculation with the STAR data [18]. We seethe results at the scale Q = 0 . p T agree well with thedata of ρ yield. In the following calculations we will fix Q = 0 . p T to provide a good p+p baseline.A hot and dense QCD matter is created shortly afterthe high energy central nucleus-nucleus collisions. Before (GeV/c) T p4 6 8 10 12 14 16 18 20 -10 -9 -8 -7 -6 -5 -4 -3 best fit) T (Q=0.5 P T Q=0.5-2.0 p wih PDF uncertainties T Q=0.5 PSTAR 200 GeV +X ρ p+p -> ) c - p ( m b G e V / d σ E d FIG. 2: Numerical calculation of the ρ production in p + pcollisions at RHIC 200 GeV comparing with STAR [18] data. a fast parton fragmented into identified hadrons in thevacuum, it should suffer energy loss due to multiple scat-tering with other partons in QCD medium. In highertwist approach, the multiple scattering is described bytwist-4 processes of hard scattering and will lead to effec-tive medium-modification of the vacuum FFs [14–16, 27–29, 32]:˜ D hq ( z h , Q ) = D hq ( z h , Q ) + α s ( Q )2 π Z Q dℓ T ℓ T × Z z h dzz h ∆ γ q → qg ( z, x, x L , ℓ T ) D hq ( z h z , Q )+ ∆ γ q → gq ( z, x, x L , ℓ T ) D hg ( z h z , Q ) i , (3)where ∆ γ q → qg ( z, x, x L , ℓ T ) and ∆ γ q → gq ( z, x, x L , ℓ T ) =∆ γ q → qg (1 − z, x, x L , ℓ T ) are the medium modified split-ting functions [27–29]. Though the medium-modifiedFFs include a contribution from gluon radiation in theQCD medium, they obey QCD evolution equations sim-ilar to the DGLAP equations for FFs in vacuum. Inthis formalism, we convolute the medium-induced kernel∆ γ q → qg ( z, x, x L , ℓ T ) and ∆ γ q → gq (instead of those vac-uum splitting functions) with the (DGLAP) evolved FFsat scale Q . We average the above medium modified frag-mentation functions over the initial production positionand jet propagation direction, scaled by the number ofbinary nucleon-nucleon collisions at the impact param-eter b in A + A collisions to replace the vacuum frag-mentation functions in Eq. (1). In the medium modifiedsplitting functions ∆ γ q → qg,gq , we can extract the depen-dency of the properties of the medium into the jet trans-port parameter ˆ q which defined as the average squaredtransverse momentum broadening per unit length. In thehigher-twist approach, the jet transport parameter ˆ q isrelated to the gluon distribution density of the medium. Phenomenologically the jet transport parameter can beassumed to be proportional to the local parton densityin the QGP phase and also to the hadron density in thehadronic gas phase [14]:ˆ q ( τ, r ) = (cid:20) ˆ q ρ QGP ( τ, r ) ρ QGP ( τ ,
0) (1 − f ) + ˆ q h ( τ, r ) f (cid:21) · p µ u µ p , (4) ρ QGP is the parton (quarks and gluon) density in an idealgas at a given temperature, f ( τ, r ) is the fraction of thehadronic phase as a function of space and time, ˆ q is thejet transport parameter at the center of the bulk mediumin the QGP phase at the initial time τ , p µ is the fourmomentum of the jet and u µ is the four flow velocity inthe collision frame.The space-time evolution of the QCD medium is givenby a full three-dimensional (3+1D) ideal hydrodynamicsdescription [33, 34]. Parton density, temperature, frac-tion of the hadronic phase and the four flow velocity atevery time-space points are provided by the hydro dy-namical model. The only free parameter is ˆ q τ , theproduct of initial value of the jet transport parameterˆ q and the time τ when the QCD medium is initiallyformed. This parameter controls the strength of jet-medium interaction, and the amount of the energy loss ofthe energetic jets. In the calculations, we use the valuesof ˆ q τ extracted in the previous studies [14–16], whichgive very nice descriptions of single π and η productionsin HIC. Moreover, we have used the EPS09 parametriza-tion sets of nuclear PDFs f a/A ( x a , µ ) to consider theinitial-state cold nuclear matter effects [35]. AA R /fm =1.2GeV @RHIC q ρ /fm =1.2GeV @RHIC q π PHENIX 0-10 % π STAR 0-10 % ρ (GeV) T p ± π AA / R ρ AA R /fm =1.2GeV qSTAR data FIG. 3: Top panel: Numerical calculation of the ρ and π production suppression factors in 0 −
10% Au+Au collisions atRHIC 200 GeV at NLO as functions of p T , comparing withSTAR [18] and PHENIX [36] data; Bottom panel: doubleratio calculation of R ρ AA /R π ± AA both at NLO, also comparingwith STAR data. (GeV) T p4 6 8 10 12 14 16 18 F r a c t i on ρ p+p 200 GeV q ρ p+p 200 GeV g ρ Au+Au 200 GeV q ρ Au+Au 200 GeV
FIG. 4: Gluon and quark contribution fraction of the totalyield both in p+p and Au+Au at RHIC
Now we are ready to calculate the single ρ productionsin heavy ion collisions up to the NLO. The nuclear mod-ification factor R AA as a function of p T is calculated todemonstrate the suppression of the production spectrumin A + A collisions relative to that in p + p collision: R AB ( b ) = dσ hAB /dyd p T N ABbin ( b ) dσ hpp /dyd p T (5)In the 0 −
10% most central Au+Au collisions at RHIC200 GeV, we calculate ρ productions at typical values ofˆ q = 1 . /fm and τ = 0 . ρ mesonat large p T region (see the top panel of Fig. 3). The the-oretical calculation and the experimental data of the π nuclear suppression factor are also presented for compar-ison. We note that the nuclear suppression factor of ρ issimilar to the one of π , as demonstrated by the doubleratio R ρ AA /R π ± AA in the bottom panel of Fig. 3, which isaround 1 calculated at the NLO accuracy. We also findthat the theoretical curve undershoots the experimentaldata of R AA same as the case in π , and the uncertaintycaused by this undershooting will be cancelled out to alarge extent when we discussing the double ratio of ρ and charged π . Here π ± FFs in vacuum are given byAKK08 [38].To understand better the nature of the suppressionpattern of ρ , we calculate the gluon (quark) contributionfraction of the total yield both in p+ p and Au+ Au colli-sions in Fig. 4. It is similar to η and π productions whichdemonstrate the domination of the quark fragmentationprocess contribution at high p T region either in p + p orin A+A collisions , and the jet quenching effect may sup-press the gluon fragmenting contribution but enhance thequark contribution. Therefore the crossing point wherethe fractional contributions of quark and gluon fragmen-tation are equal, will move toward lower p T in Au + Aucollision, as one observes in Fig. 4.
20 40 60 80 -11 -10 -9 -8 -7 -6 -5 -4 +X 2.76 TeV ρ Pb+Pb->=0.6 fm τ ρ Pb+Pb->=0.6 fm τ /fm ± =2.2 q ρ
20 40 60 80
CMS Charged hadron 0-10%/fm ± =2.2 q ρ /fm ± =2.2 q π (GeV) T p ) c - p ( m b G e V / d σ E d AA R FIG. 5: Numerical calculation of the ρ production in 0 −
10% Pb + Pb collisions at LHC 2 .
76 GeV in the top panel;theoretical calculation results of nuclear suppression factor of ρ and π are compared with the experimental data of chargedhadron [37] in 0 −
10% Pb + Pb collisions at LHC 2 .
76 GeVin the bottom panel.
We also predict the ρ production in the 0 −
10% mostcentral Pb + Pb collisions at the LHC with √ s NN =2 .
76 TeV in the top panel of Fig. 5. The values of the ˆ q are set to be the same as the typical values which havebeen used to describe production suppression of both sin-gle π and η mesons at the LHC [14–16]. We can see that,with the increase of p T , the nuclear modification factor of ρ meson goes up slowly. In the calculation, best fit to thePHENIX data on π nuclear suppression factor as a func-tion of p T in 0 −
5% Au+Au collisions at √ s = 200 GeVgives ˆ q = 1 . ± .
30 GeV / fm. Similarly, the best fitto the CMS data on charged hadron nuclear suppressionfactor in 0 −
5% Pb+Pb collisions at √ s = 2 .
76 TeV asa function of p T would gives ˆ q = 2 . ± . / fm at τ = 0 . / c [13]. The same values of ˆ q τ have beenemployed to give a very nice description of ρ productionsin LHC shown in the bottom panel of Fig. 5.To compare the different trends of π and ρ spectra,we plot the ratio ρ /π as a function of the transversemomentum p T in Fig. 6. As we have mentioned that inthe study the π FFs are given by AKK08 [38]. We notethat even the validity of the π (charged hadron) FFshad been challenged by the over-predicting of its pro-duction in the LHC and Tevatron due ot the too-hardgluon-to-hadron FFs in the parameterizations [39]. A re-cent attempt to address the problem and a global refitis performed in Ref. [40]. The uncertainty brought in bythe usage of AKK08 fortunately do not affect the resultsof the nuclear modification factor R AA much due to thecancellation when taking the ratio of A+A production top+p reference. Therefore one expects that the extrac-tion of jet transport parameter ˆ q from the comparisonbetween theoretical calculated R AA and the experimentaldata will not be affected much by such FFs uncertainties.In the studies of particle ratio, π fragmentation functionand its jet chemistry are used as reference to understandother mesons such as η , its FFs uncertainties certainlywill be expected to affect particle ratios like η/π . How-ever, since light mesons such as π and η are dominatedby quark fragmenting contribution, such effect is there-fore minimized.Fig. 6 illustrates that the ratio ρ /π increases withthe p T in p + p collision at the RHIC energy and LHC.Though the jet quenching effect may alter the ratio a lit-tle bit in A + A at lower p T , as p T becomes larger, theratio in A+A comes very close that in p+p, especially atthe LHC with higher p T . We note flat curves are observedin η/π ratios as functions of p T at both the RHIC andthe LHC, whereas a increasing ρ /π with respect to p T are shown in Fig. 6. The ρ /π ratio in the RHIC demon-strates a more rapidly increasing behavior with respectto p T . It is realized that the flat particle ratio depen-dence of p T is therefore not a universal trend, and theshape of the particle ratio depends on the relative slopeof their spectra in p+p , different flavor contributions toFFs as well as flavor dependence of parton energy loss inthe QGP. (GeV) T p10 20 30 40 50 60 70 80 90 r a t i o π / ρ p+p @LHC π / ρ Au+Au @RHIC π / ρ Pb+Pb @LHC π / ρ FIG. 6: ρ /π production ratio as a function of final state p T calculated both in p+p and A+A collisions at RHIC andLHC We note that at high p T region, the productions ofboth ρ and π are dominated by quark contribution (forexample, see Fig. 4). If at high p T , quark FFs of ρ and π have a relatively weak dependence on z h and p T , then we have:Ratio( ρ /π ) = dσ η dp T / dσ π dp T ≈ R F q ( p T z h ) D q → ρ ( z h , p T ) dz h z h R F q ( p T z h ) D q → π ( z h , p T ) dz h z h ≈ Σ q D q → ρ ( h z h i , p T )Σ q D q → π ( h z h i , p T ) . Therefore, while quark and gluon may lose different frac-tions of their energies, at very high p T region, the ratio ρ /π in A + A collisions should approximately be de-termined only by quark FFs in vacuum with the p T shiftbecause of the parton energy loss. As we can see in Fig. 1,the quark FFs at large scale Q (= p T ) change slowly withthe variation of both z h and p T , then the ratio of ρ /π in both A + A and p + p may approach to each other atlarger p T . It is just as we have observed in the case forthe yield ratio of η/π [16]. Acknowledgments:
This research is supported bythe MOST in China under Project No. 2014CB845404,NSFC of China with Project Nos. 11435004, 11322546,11521064, and partly supported by the Fundamental Re-search Funds for the Central Universities, China Univer-sity of Geosciences (Wuhan) (No. 162301182691) [1] X. N. Wang and M. Gyulassy, Phys. Rev. Lett. , 1480(1992).[2] M. Gyulassy, I. Vitev, X. N. Wang and B. W. Zhang, In*Hwa, R.C. (ed.) et al.: Quark gluon plasma* 123-191[nucl-th/0302077].[3] K. Aamodt et al. [ALICE Collaboration], Phys. Rev.Lett. , 092301 (2012) [arXiv:1110.0121 [nucl-ex]].[4] C. Adler et al. [STAR Collaboration], Phys. Rev. Lett. , 082302 (2003) [nucl-ex/0210033].[5] A. Adare et al. [PHENIX Collaboration], Phys. Rev. C , 024908 (2009) [arXiv:0903.3399 [nucl-ex]].[6] B. I. Abelev et al. [STAR Collaboration], Phys. Rev. C , 034909 (2010) [arXiv:0912.1871 [nucl-ex]].[7] I. Vitev, S. Wicks and B. W. Zhang, JHEP , 093(2008) [arXiv:0810.2807 [hep-ph]].[8] I. Vitev and B. W. Zhang, Phys. Rev. Lett. , 132001(2010) [arXiv:0910.1090 [hep-ph]].[9] W. Dai, I. Vitev and B. W. Zhang, Phys. Rev. Lett. ,no. 14, 142001 (2013) [arXiv:1207.5177 [hep-ph]].[10] G. Aad et al. [ATLAS Collaboration], Phys. Rev. Lett. , 252303 (2010) [arXiv:1011.6182 [hep-ex]].[11] S. Chatrchyan et al. [CMS Collaboration], Phys. Rev. C , 024906 (2011) [arXiv:1102.1957 [nucl-ex]].[12] Z. B. Kang, R. Lashof-Regas, G. Ovanesyan, P. Saadand I. Vitev, Phys. Rev. Lett. , no. 9, 092002 (2015)[arXiv:1405.2612 [hep-ph]].[13] K. M. Burke et al. [JET Collaboration], Phys. Rev.C , no. 1, 014909 (2014) [arXiv:1312.5003 [nucl-th]];Z. Q. Liu, H. Zhang, B. W. Zhang and E. Wang, Eur.Phys. J. C , no. 1, 20 (2016) [arXiv:1506.02840 [nucl-th]].[14] X. F. Chen, C. Greiner, E. Wang, X. N. Wang and Z. Xu,Phys. Rev. C , 064908 (2010) [arXiv:1002.1165 [nucl-th]].[15] X. F. Chen, T. Hirano, E. Wang, X. N. Wangand H. Zhang, Phys. Rev. C , 034902 (2011) [arXiv:1102.5614 [nucl-th]].[16] W. Dai, X. F. Chen, B. W. Zhang and E. Wang, Phys.Lett. B , 390 (2015) [arXiv:1506.00838 [nucl-th]].[17] W. Dai and B. W. Zhang, arXiv:1612.05848 [hep-ph].[18] G. Agakishiev et al. [STAR Collaboration], Phys. Rev.Lett. , 072302 (2012) [arXiv:1110.0579 [nucl-ex]].[19] A. Adare et al. [PHENIX Collaboration], Phys. Rev. C , 024909 (2011) [arXiv:1004.3532 [nucl-ex]].[20] R. Bala, I. Bautista, J. Bielcikova and A. Ortiz,Int. J. Mod. Phys. E , no. 07, 1642006 (2016)[arXiv:1605.03939 [hep-ex]].[21] W. Liu, C. M. Ko and B. W. Zhang, Phys. Rev. C ,051901 (2007) [nucl-th/0607047].[22] S. J. Brodsky and A. Sickles, Phys. Lett. B , 111(2008) [arXiv:0804.4608 [hep-ph]].[23] X. Chen, H. Zhang, B. W. Zhang and E. Wang, J. Phys. , 015004 (2010) [arXiv:0806.0556 [hep-ph]].[24] H. Saveetha, D. Indumathi and S. Mitra, Int. J. Mod.Phys. A , no. 07, 1450049 (2014) [arXiv:1309.2134[hep-ph]].[25] D. Indumathi and H. Saveetha, Int. J. Mod. Phys. A ,1250103 (2012) [arXiv:1102.5594 [hep-ph]].[26] M. Hirai and S. Kumano, Comput. Phys. Commun. ,1002 (2012) [arXiv:1106.1553 [hep-ph]].[27] X. f. Guo and X. N. Wang, Phys. Rev. Lett. , 3591(2000) [hep-ph/0005044].[28] B. W. Zhang and X. N. Wang, Nucl. Phys. A , 429(2003) [arXiv:hep-ph/0301195].[29] B. W. Zhang, E. Wang and X. N. Wang, Phys. Rev.Lett. , 072301 (2004) [nucl-th/0309040]; A. Schafer,X. N. Wang and B. W. Zhang, Nucl. Phys. A , 128 (2007) [arXiv:0704.0106 [hep-ph]].[30] N. Kidonakis and J. F. Owens, Phys. Rev. D , 054019 (2001) doi:10.1103/PhysRevD.63.054019[hep-ph/0007268].[31] H. L. Lai et al. [CTEQ Collaboration], Eur. Phys. J. C , 375 (2000) [hep-ph/9903282].[32] W. Dai, B. W. Zhang, H. Z. Zhang, E. Wang andX. F. Chen, Eur. Phys. J. C , no. 8, 571 (2017)[arXiv:1702.01614 [nucl-th]].[33] T. Hirano, Phys. Rev. C , 011901 (2002).[34] T. Hirano and K. Tsuda, Phys. Rev. C , 054905 (2002).[35] K. J. Eskola, H. Paukkunen and C. A. Salgado, JHEP , 065 (2009) [arXiv:0902.4154 [hep-ph]].[36] S. S. Adler et al. [PHENIX Collabora-tion], Phys. Rev. Lett. , 072301 (2003)doi:10.1103/PhysRevLett.91.072301 [nucl-ex/0304022].[37] K. Aamodt et al. [ALICE Collaboration], Phys. Lett.B , 30 (2011) doi:10.1016/j.physletb.2010.12.020[arXiv:1012.1004 [nucl-ex]].[38] S. Albino, B. A. Kniehl and G. Kramer, Nucl. Phys.B , 42 (2008) doi:10.1016/j.nuclphysb.2008.05.017[arXiv:0803.2768 [hep-ph]].[39] D. d’Enterria, K. J. Eskola, I. Helenius andH. Paukkunen, Nucl. Phys. B , 615 (2014)doi:10.1016/j.nuclphysb.2014.04.006 [arXiv:1311.1415[hep-ph]].[40] D. de Florian, R. Sassot, M. Epele, R. J. Hernndez-Pintoand M. Stratmann, Phys. Rev. D91