Measurements of the transverse-momentum-dependent cross sections of $J/\psi$ production at mid-rapidity in proton+proton collisions at $\sqrt{s} =$ 510 and 500 GeV with the STAR detector

We present measurements of the differential production cross sections of the inclusive $J/\psi$ meson as a function of transverse momentum ($p_{T}^{J/\psi}$) using the $\mu^{+}\mu^{-}$ and $e^{+}e^{-}$ decay channels in proton+proton collisions at center-of-mass energies of 510 and 500 GeV, respectively, recorded by the STAR detector at the Relativistic Heavy Ion Collider. The measurement from the $\mu^{+}\mu^{-}$ channel is for 0 $<p_{T}^{J/\psi}<$ 9 GeV/$c$ and rapidity range $|y^{J/\psi}|<$ 0.4, and that from the $e^{+}e^{-}$ channel is for 4 $<p_{T}^{J/\psi}<$ 20 GeV/$c$ and $|y^{J/\psi}|<$ 1.0. The $\psi(2S)$ to $J/\psi$ ratio is also measured for 4 $<p_{T}^{\rm meson}<$ 12 GeV/$c$ through the $e^{+}e^{-}$ decay channel. Model calculations, which incorporate different approaches toward the $J/\psi$ production mechanism, are compared with experimental results and show reasonable agreement within uncertainties. A more discriminating comparison to theoretical models at low $p_T$ can be performed in the future, if the calculations are carried out within our fiducial volume, eliminating the uncertainty due to the $J/\psi$ polarization.


I. INTRODUCTION
The J/ψ meson is a bound state of charm and anticharm quarks (cc) which was discovered several decades ago [1]. In hadronic collisions at energies reached at the Relativistic Heavy Ion Collider (RHIC), J/ψ are primarily produced via inelastic scattering by two gluons into charm and anti-charm quarks, followed by hadronization of the cc pair [2,3]. Studying the J/ψ production provides valuable knowledge for the understanding of Quantum Chromodynamics (QCD) in both perturbative and non-perturbative regimes. The production of the cc pair can be calculated using the perturbative approach, however, the evolution of a cc pair into a J/ψ meson is nonperturbative, and the theoretical description remains a challenge. Different theoretical approaches have been proposed to describe the J/ψ production mechanism [4][5][6][7][8]. However, these descriptions have difficulties in explaining the experimental results of production cross section and polarization simultaneously. Therefore, precise measurements of the J/ψ cross section in elementary collisions at different collision energies are essential for investigating the J/ψ production mechanism. Moreover, as an important probe of the hot and dense medium, known as quark-gluon plasma (QGP), it is necessary to have a good understanding of the J/ψ production mechanism in elementary collisions in order to help understand the modification to its production in heavy-ion collisions, which has been proposed and widely pursued to study the properties of QGP [9].
There are three notable models for J/ψ production which differ mainly in the description of the non-perturbative process. These are the Color Singlet Model (CSM) [4], Non-Relativistic QCD formalism (NRQCD) [5] and the Color Evaporation Model (CEM) [6]. In the CSM, it is assumed that the hadronization process does not change the quantum numbers of the cc pair. The initially produced cc can then bind to a given charmonium state only if it is created in a color-singlet state with matching angular-momentum quantum numbers. The Next-to-Next-to-Leading-Order CSM (NNLO ⋆ CSM) has been tested for the S-wave quarkonium states in the Tevatron and LHC data. However, this model is not able to calculate the full NNLO contribution, or provide the predictions for the P -wave states due to techni-cal issues [10]. Therefore, it is expected to underestimate the production for the quarkonium states which have significant contributions from the decays of excited states, known as feed-down contributions [11]. In the NRQCD approach, the charmonium can be produced from both the CS state and a color-octet (CO) state. The color neutralization of the CO state is achieved by radiating soft gluons during the hadronization process. In the CEM, the produced cc pair is assumed to evolve into a J/ψ with a certain probability if its invariant mass is below the threshold for producing a DD pair. In this model, spin is always summed over which prevents it from predicting the J/ψ polarization. A recent improvement to the CEM (Improved CEM (ICEM) [7]) overcomes this issue by sorting out different spin states and is able to predict the polarization of the quarkonium states. In the low transverse momentum (p T ) range of the charmonium, the cc cross section becomes difficult to calculate at collider center-of-mass energies since the dynamics are sensitive to the large logarithms of small Bjoken x. A newly developed Color Glass Condensate (CGC) effective theory of small-x QCD provides a viable path towards calculating the J/ψ cross section at low p T (p T < ∼M , where M is the quarkonium mass) by combining the CGC effective theory with the NRQCD formalism [8].
This paper presents the measurements of the J/ψ production cross sections covering a wide p T range from 0 to 20 GeV/c in proton+proton collisions at center-ofmass energies of 510 and 500 GeV at RHIC. These cross sections are measured in two decay channels, which include the µ + µ − channel, for 0 < p J/ψ T < 9 GeV/c and J/ψ rapidity (|y J/ψ |) < 0.4, and the e + e − channel, for 4 < p J/ψ T < 20 GeV/c and |y J/ψ | < 1.0, respectively. The measured cross sections contain the direct production of J/ψ, contributions from excited charmonium states, and from decays of bottom-flavored hadrons. The first two are often categorized as prompt J/ψ as they are produced at the collision vertex and cannot be experimentally separated. The last one is often called non-prompt J/ψ, while the detector setup used in this analysis cannot experimentally distinguish it from prompt J/ψ. The feed-down contribution to the J/ψ production is an additional complication in understanding the J/ψ production mechanism as nearly 30 − 40% of the inclusive J/ψ yields come from the decay of excited charmonium states [12,13].
Many experiments have already presented the results of heavy quarkonium production in electron+positron, hadron+hadron and heavy-ion collisions [14,15]. The latest measurements from the LHC [16][17][18] probe the high p T production cross sections in proton+proton collisions with energies of 7, 8, and 13 TeV. The large kinematic range of the J/ψ measurement at the highest beam energies at RHIC (510 and 500 GeV) provides valuable insights to the J/ψ production mechanism. Additionally, the ψ(2S) to J/ψ ratio is measured in the e + e − decay channel in the p T range of 4 < p meson T < 12 GeV/c. This measurement could help constrain the feed-down contribution to the J/ψ from the excited charmonium states.
The paper is organized as follows: the STAR detector will be discussed in section II, and the analyses of the µ + µ − and e + e − decay channels will be described in detail in section III and IV, respectively. The results from these two different channels will be presented and compared to different theoretical models in section V. Finally, conclusions will be given in section VI.

II. THE STAR DETECTOR
The STAR detector is optimized for high energy nuclear physics. It has excellent particle identification capability and a large acceptance at mid-rapidity. The heart of the STAR detector is the Time Projection Chamber (TPC). The TPC is the primary tracking detector for charged particles and provides particle identification via measurements of their ionization energy losses (dE/dx) [19]. It covers the full azimuthal range (0 ≤ φ < 2π) and a large pseudorapidity range (|η| < 1). The p T of charged particles are measured from the curvature of their trajectories in the 0.5 Tesla solenoidal field [20]. There are 30 iron bars, known as "backlegs" outside the coil to provide the return flux path for the magnetic field. These are 61 cm thick at a radius of 363 cm corresponding to about 5 interaction lengths. These backlegs play an essential role in enhancing the muon purity by absorbing the background hadrons from collisions. The hadron rejection rate is about 99% as shown in the simulation study [21]. The Muon Telescope Detector (MTD) is a fast detector which uses Multi-gap Resistive Plate Chamber technology to record signals, also referred to as "hits", generated by charged particles traversing it. It provides single-muon and dimuon triggers depending on the number of hits within a predefined online timing window. The MTD modules are installed at a radius of about 403 cm, and the full MTD detector covers about 45% in azimuth within |η| < 0.5 [21]. The timing resolution of the MTD is ∼100 ps and the spatial resolutions are ∼1-2 cm in both rφ and z directions as demonstrated in cosmic-ray data [22]. The data used in this analysis were taken during the run in which the MTD detector was 63% completed. The Barrel Electromagnetic Calorime-ter (BEMC) is a lead-scintillator sampling calorimeter with 23 radiation lengths [23]. The BEMC, being a thick absorber, is dedicated to measuring energies of particles with electromagnetic interactions, such as electrons and positrons. The BEMC is physically segmented into a total of 4800 towers with a granularity of 0.05 × 0.05 in ∆φ×∆η. The energy deposited in the towers is used as a trigger to record rare events. The Vertex Position Detectors (VPD) [24] and the Beam Beam Counters (BBC) [25] are scintillator-based detectors located on both sides of the main detector, and they cover pseudorapidity ranges from 4.4 to 4.9 and 2 to 5, respectively.

A. Data and Monte Carlo
Data for the µ + µ − channel in this analysis were collected by the STAR detector during the 2013 RHIC pro-ton+proton run at a collision energy of 510 GeV. The corresponding integrated luminosity sampled by the MTD dimuon trigger, which requires at least two coincidence hits on the MTD, as well as signals in the VPD and the two BBCs, within the bunch crossing is 22.0 pb −1 . Events used in the analysis are required to have a valid reconstructed vertex with at least two tracks that are associated with corresponding MTD hits.
A Monte Carlo (MC) simulation sample was generated by a single-particle generator with flat distributions in p T , φ and y for the J/ψ signal. These simulated signals were passed through a full GEANT3 [26] STAR detector simulation and then "embedded" into real data events. These embedded events were reconstructed using the same reconstruction procedure used for real data. The kinematic distributions of the embedded J/ψ were weighted by the p T spectrum of J/ψ in proton+proton collisions at a collision energy of 510 GeV, determined via interpolation through a global fit of world-wide measurements of J/ψ cross sections [27]. Due to the systematic uncertainties on various distortion corrections for the TPC, the p T resolution of the reconstructed muon in MC was retuned to match the reconstructed J/ψ signal mass shape in data.

B. Muon candidate selection
Muon candidates for reconstructing the J/ψ signal must satisfy the following selection criteria: p µ T is greater than 1.3 GeV/c to ensure the track can reach the MTD detector; the pseudorapidity of the track is within the MTD acceptance, |η µ | < 0.5; the distance of closest approach (DCA) to the collision vertex must be less than 3 cm to suppress background tracks from pile-up events and secondary decay vertices; the number of TPC clusters used in track reconstruction is more than 15 (the maximum possible is 45) to ensure good momentum resolution; the number of TPC clusters used for the dE/dx measurement should be more than 10 to have good dE/dx resolution; the ratio of the number of TPC clusters used over the number of possible clusters is at least 0.52 to avoid double counting for the same tracks from track splitting. Tracks are propagated from the interaction vertex to the MTD and required to match the MTD hits geometrically which fired the trigger. In addition, the muon candidates were selected by an advanced muon identification method called the Likelihood Ratio method which is described in Ref. [28]. The rapidity of the µ + µ − pairs should be smaller than 0.4 to reduce the edge effect from the J/ψ kinematic acceptance which will be described in the next section. Figure 1 shows the invariant mass spectrum of the µ + µ − pairs with the selection criteria described above applied to both candidate daughters. This can be well described by a single Gaussian as signal plus second-order polynomial function as background. A total of 1154 ± 54 final J/ψ candidates are observed within the kinematic phase space of 0 < p J/ψ T < 9 GeV/c and |y J/ψ | < 0.4. C. Acceptance and efficiency The differential production cross section multiplied by the J/ψ → µ + µ − branching ratio (BR), (5.961 ± 0.033)% [29], is given by where N corr. J/ψ→µ + µ − is the efficiency-corrected number of J/ψ candidates. Ldt is the corresponding integrated luminosity. ∆p T and ∆y are the corresponding bin widths in p T and y of the µ + µ − pairs, respectively. For each µ + µ − candidate, a weighting factor (w i ) is multiplied to correct for the detector acceptance (A) and the total reconstruction efficiency (ε reco. ), to obtain N corr.
and i indicates the ith candidate. This minimizes the potential bias from the efficiency correction due to the gaps between the MTD modules and restricted pseudorapidity coverage.
The detector acceptance, A, is the probability of detecting muons having certain kinematics, namely p µ T > 1.3 GeV/c and |η| < 0.5, from the J/ψ decay within the detector fiducial volume. The A can be factorized into the J/ψ decay kinematic acceptance and the MTD geometric acceptance, A = A J/ψ × A MTD . The acceptance can be determined with the muon angular distribution calculated in the J/ψ rest frame by the following formula [30]: where θ ⋆ is the polar angle between the µ + momentum in the J/ψ rest frame and the direction of the J/ψ momentum in the laboratory frame; φ ⋆ is the azimuthal angle between the J/ψ production plane (defined in the J/ψ rest frame by the momenta of the incoming protons) and the J/ψ decay plane in the lab frame; and λ i are the parameters for different polarization configurations. Similar to the analyses carried out by other experiments [11,16], five extreme configurations are considered to cover the polarization phase space: unpolarized, λ θ = λ φ = λ θφ = 0; longitudinaly polarized, λ θ = −1, λ φ = λ θφ = 0; zero transversely polarized, λ θ = +1, λ φ = λ θφ = 0; positively transversely polarized, λ θ = +1, λ φ = +1, λ θφ = 0; and negatively transversely polarized, λ θ = +1, λ φ = −1, λ θφ = 0. Figure 2 shows the J/ψ decay kinematics acceptance and the MTD geometric acceptance as a function of J/ψ p T with different polarization assumptions, respectively. There is a significant difference in the J/ψ decay kinematic acceptance for different polarization assumptions at p T around 2 GeV/c, and the fractional difference becomes smaller at higher p T . On the other hand, the MTD geometric acceptance is almost independent of the J/ψ polarization configuration. The total J/ψ reconstruction efficiency, ε reco. , includes the VPD requirement, the TPC tracking efficiency, the  vertex-finding efficiency, the dimuon triggering, the MTD in-situ response and matching, and muon-identification efficiencies, as shown in the following: where the superscripts, 1 and 2, indicate the first and second muon from a J/ψ candidate. The VPD efficiency (ε VPD ) is obtained from the zero-bias MC sample. The TPC tracking efficiency (ε trk ) is calculated from the J/ψ → µ + µ − MC sample. The vertex finding efficiency (ε vtx ) is obtained from data directly and is about 95% across the entire J/ψ p T region. The MTD trigger efficiency (ε trig. ) includes the trigger electronics efficiency which varies from 95% at low p T to more than 99% at high p T , and the online timing window cut efficiency which reaches a plateau of 99.9%. The MTD efficiency (ε MTD ) is determined from the cosmic ray data for the in-situ response efficiency and from the MC sample for the matching efficiency. It is evaluated as a function of muon p T for each MTD backleg and module separately. Finally, the muon identification efficiency is calculated from the MC events and the plateau efficiency is above 95% [28]. Figure 3 shows the individual and total efficiencies used for the J/ψ cross section measurement as a function of p The N corr. J/ψ→µ + µ − in different p T regions are extracted using the χ 2 fit with several combinations of signal and background models of the efficiency-corrected µ + µ − mass distributions. The signal shape is modeled by a single Gaussian, double Gaussian or Crystal-Ball function, and the background shape can be well described by the samesign muon track pairs or a polynomial function at different orders. The averaged result from the various fits with different shapes for signal and background is used as the mean of N corr.
J/ψ , and the maximum deviation from the mean is assigned as the signal extraction systematic uncertainty. Figure 4 shows an example of the fit results using a single Gaussian as the signal function and samesign muon track pairs as the background template for different p T bins.

D. Uncertainties
The statistical uncertainty is about 10.6−18.8% for different J/ψ p T regions. There are several systematic uncertainties considered in this analysis. The maximum deviation in the results for variations in cuts/methods is taken as the systematic uncertainty for each source listed below.
• The luminosity estimate and the in-bunch pile-up effect contribute global uncertainties of 8.1% [31] and 7.7%, respectively. The latter is estimated by comparing the results with different selections of the difference of the z position measured by the VPD and TPC detectors.
• The uncertainty in the number of J/ψ candidates extraction is evaluated using different signal and background models as described in the previous section. It contributes about 0.3−4.8% uncertainty depending on J/ψ p T .
• The uncertainty in the TPC tracking efficiency is estimated by comparing the results using different TPC track quality selection criteria and contributes an overall 10.6% uncertainty.
• The uncertainty in the MTD trigger includes two components: (i) the trigger electronics which is dominated by the statistical uncertainty in calculating the MTD trigger efficiency; (ii) the online timing window cut which is evaluated by the difference between the 2013 and 2015 data-taking. A total 3.6% uncertainty is assigned.
• The MTD efficiency includes three sources: (i) statistical precision of the cosmic ray data; (ii) different fit templates used for determining the response efficiency which is the main contributor; (iii) difference in the matching efficiency between cosmic ray data and simulation. The resulting uncertainty is between 1.9−7.6%.
• The uncertainty in the muon identification efficiency is determined by comparing the efficiencies from the data-driven method and the MC sample. It contributes a 5.2−8.7% uncertainty depending on J/ψ p T [28].
• The uncertainty in vertex finding efficiency is estimated by comparing the efficiencies from datadriven and zero-bias MC sample. A 4.1% uncertainty is assigned.
The systematic uncertainties from the MTD geometric acceptance is negligible. Systematic uncertainties from different sources are added in quadrature. Figure 5 shows all the uncertainties as a function of J/ψ p T .
E. Cross-section for J/ψ → µ + µ − We measure the invariant differential cross section multiplied by the µ + µ − branching ratio of the J/ψ meson as a function of J/ψ p T in a fiducial volume defined by p µ T > 1.3 GeV/c and |η µ | < 0.5 (fiducial cross section) and in a full muon decay phase space with |y J/ψ | < 0.4 (full cross section). Figure 6 shows the fiducial and full cross sections of the J/ψ production. The fiducial cross section is calculated using the fiducial weight, w fid. The proton+proton collision data at √ s = 500 GeV, used in the e + e − analysis, were recorded by the STAR detector in 2011. The integrated luminosity of the data set is 22.1 pb −1 sampled by the BEMC trigger which requires a BEMC tower with a transverse energy deposit larger than 4.3 GeV [14]. The e ± candidates are reconstructed and identified using information from the TPC and BEMC detectors. The track quality requirements are that each track has at least 25 out of 45 possible hits in the TPC, the number of hits for the dE/dx measurement must be larger than 15, and tracks are reconstructed within the TPC acceptance of |η| < 1. The electron and positron candidates are then identified by their ionization energy loss ( dE/dx ) in the TPC. More than 15 TPC hits are required to calculate dE/dx to ensure a good dE/dx resolution. The normalized dE/dx is  2πp T dp T dy ± δstat. ± δsys. BR× dσ 2 full 2πp T dp T dy ± δstat. ± δsys. 5.68 (6.7 ± 0.9 ± 0.9) ×10 0 (2.6 ± 0.3 ± 0.3 +4.3 −0.6 ) ×10 1 7.0 -9.0 7.73 (1.7 ± 0.3 ± 0.2) ×10 0 (4.6 ± 0.9 ± 0.6 +3.7 −0.9 ) ×10 0 TABLE I. A summary of the fiducial and full differential cross sections of the J/ψ production via the µ + µ − decay channel in proton+proton collisions at √ s = 510 GeV. The fiducial volume is defined as p µ T > 1.3 GeV/c and |η µ | < 0.5. A common luminosity uncertainty of 11.2% is not included. defined as follows: where dE/dx me. and dE/dx exp. e are the measured dE/dx and the expected dE/dx value for electron, and the σ dE/dx is the experimental ln(dE/dx) resolution. The nσ e requirement for the triggered e ± candidates is set to be |nσ e | < 2. The triggered e ± candidate is also required to have p T > 3.5 GeV/c, its track must have DCA from the primary vertex less than 1 cm, and it is required to match to a BEMC trigger tower in which the ADC value is larger than 290, corresponding to a deposited energy of 4.3 GeV. A cut on the ratio of the momentum measured by the TPC to the energy deposited in the BEMC towers, 0.3 < pc/E < 1.5, is used to further suppress contribution of hadrons in triggered electron selection. However, non-triggered e ± candidates only require TPC tracks with additional requirements of p T > 1 GeV/c and DCA < 3 cm.   Figure 7 shows both J/ψ and ψ(2S) signals which are reconstructed from e + e − pairs, which include the triggered electron pairing with either another triggered electron or with a non-triggered electron. The signal shape for J/ψ candidates is obtained from the MC simulation, which includes track momentum resolution and electron bremsstrahlung radiation in the detector. On the other hand, the background includes the combinatorial background evaluated using the like-sign e + e + and e − e − pairs within the same event (green histogram) and the residual background, which mainly comes from the Drell-Yan process, cc and bb decays. The invariant mass distribution after the like-sign background subtraction is fitted with a J/ψ signal shape combined with an exponential function. The raw J/ψ yield is obtained by counting the bin contents after subtracting the residual background in the mass range of 2.7 < M ee < 3.3 GeV/c 2 . There are 9581 ± 207 J/ψ signals in 4 < p J/ψ T < 20 GeV/c. About ∼10% of J/ψ candidates are reconstructed outside this mass window based on MC simulations, and the J/ψ raw yield as a function of p T is corrected for this effect. A total of 350 ± 89 ψ(2S) signals are obtained in the mass counting range of 3.5 < M ee < 3.8 GeV/c 2 . The ψ(2S) signals are extracted by using the same method as for J/ψ, but with a linear function to describe the residual background. Figure 8 shows invariant mass distributions of e + e − pairs (black histogram) and the like-sign (filled histogram) e + e + and e − e − pairs pairs in three represen-tative p T bins.
The measurement of the J/ψ differential production cross section multiplied by BR for the e + e − decay channel, (5.971 ± 0.032)% [29], is defined as: where BR is the branching ratio for the J/ψ → e + e − decay channel; N raw J/ψ→e + e − is the raw number of reconstructed J/ψ via the e + e − pairs; Aε is the detector's geometric acceptance times the detection efficiency of the J/ψ candidates; Ldt is the corresponding integrated luminosity; and ∆p T and ∆y are the corresponding bin widths in p T and rapidity of the e + e − pairs, respectively.
The total J/ψ detection efficiency, Aε, includes the detector acceptance, the mass bin counting efficiency, and individual efficiency of the electron candidates including the TPC tracking efficiency, the electron identification, the selection on the number of hits for the dE/dx measurement, the additional efficiency for the trigger, and the cut on pc/E for the triggered electrons. The decay electron's momentum resolution and additional p T smearing are also included in the calculation of the J/ψ detection efficiency. The efficiencies for the number of dE/dx hits and nσ e cuts are assessed using a pure electron sample from photon conversion data, while all the other acceptance and efficiencies are obtained from MC simulation with the STAR detector geometry. Figure 9 shows the individual efficiencies for the triggered electron candidates as a function of p e T . The detection efficiency for ψ(2S) candidates is obtained in the same way as for J/ψ. The relatively larger invariant mass of ψ(2S) enhances the trigger efficiency in the low p T range while the slightly larger opening angle between the electron and positron daughters decaying from the ψ(2S) will result in a smaller acceptance and thus a lower detection efficiency at high p T . Figure 10 shows the detection efficiency for J/ψ and ψ(2S), as well as the ψ(2S) to J/ψ efficiency ratio as a function of p e + e − T .

B. Uncertainties
The systematic uncertainties for the final J/ψ cross section are estimated by varying analysis selections in both data and MC simulation and comparing the corresponding J/ψ cross section to the nominal value. The systematic uncertainties considered in this analysis are the following: • The uncertainty in the luminosity contributes an overall 8.1% uncertainty [31].
• The J/ψ extraction uncertainty is estimated by using different fitting ranges and different residual background shapes. It contributes a 0.2−12.7% uncertainty depending on J/ψ p T .
• The uncertainty in the TPC tracking efficiency is estimated by comparing the results using different TPC track quality selections, and it contributes a 4−14% uncertainty.
• The uncertainty in the trigger efficiency is evaluated by comparing BEMC response in data and MC simulation, and the contribution is a 0.3−11.8% in various J/ψ p T range.
• The uncertainty in the electron identification is estimated by comparing the difference between photonic electron's nσ e distribution at different p T ranges, and it contributes an overall 1% uncertainty.
• The J/ψ internal conversion is estimated to be 4% within the mass counting range of 2.7 to 3.3 GeV/c 2 .
The systematic uncertainty from the vertex finding is negligible. Systematic uncertainties from different sources are added in quadrature. Figure 11 shows all the uncertainties as a function of J/ψ p T . For the ψ(2S) to J/ψ ratio, most of the systematic uncertainties cancel in the ratio, except for the signal extraction and trigger efficiency. They are evaluated the same way as for the J/ψ. They contribute 5.5% and 5.3%, respectively, to the uncertainty of the ratio measurement.
C. Cross-section for J/ψ → e + e − Figure 12 shows the fiducial and full production cross sections measured in the e + e − decay channel for 4 < p J/ψ T < 20 GeV/c and |y| < 1.0. Table II summarizes the cross sections of the J/ψ production as a function of p T . Similarly as for the µ + µ − channel, the full cross section is calculated under the unpolarized assumption and the polarization envelope is obtained from the same five extreme cases. This provides us the information of the J/ψ production cross section within the full J/ψ decay phase space. On the other hand, the fiducial cross section only accesses the restricted phase space, but it is independent from any polarization assumptions. The integrated fiducial and full cross sections from 4 to 20 GeV/c of J/ψ p T are 2.90 ± 0.08 (stat.) ± 0.22 (sys.) ± 0.24 (lumi.) nb and 10.7 ± 0.5 (stat.) ± 0.8 (sys.) +5.9 −2.2 (pol.) ± 0.9 (lumi.) nb, respectively.
The cross section ratio of ψ(2S) over J/ψ is 0.038 ± 0.010 (stat.) ± 0.003 (sys.) measured in the p T range of 4 < p meson T < 12 GeV/c, which is shown in Fig. 13. The result is consistent with other experimental measurements [32][33][34][35][36] and there is no obvious dependence on 2πp T dp T dy ± δstat. ± δsys. BR× dσ 2 full 2πp T dp T dy ± δstat. ± δsys.  the collision energy observed. The ICEM model [7] prediction is consistent with our data within uncertainties.    Figure 14 shows the differential cross section of inclusive J/ψ production in proton+proton collisions at √ s = 510 and 500 GeV measured by the STAR experiment combining the µ + µ − and e + e − decay channels. Please note that there is a ∼3% difference between the cross sections at 510 and 500 GeV collision energies and the difference between the different rapidity coverages is negligible [37]. The unpolarized-assumption result is compared to the NRQCD [5] and ICEM [7] calculations of J/ψ production, which includes feed-down contributions from excited charmonium states. The prediction from the CGC effective theory coupled with NRQCD (CGC+NRQCD) [8] lies systematically above the data at low p T , however, it is consistent with the data within the polarization envelope. The NLO NRQCD [5] result does a reasonably good job in describing the data above 6 GeV/c. The ICEM calculation [7] can cover the entire p T range and is also consistent with the data within the polarization envelope. The feed-down contributions from B-hadrons are about 10−20% in the p T < 10 GeV/c region and nearly 40% in our maximum p T bin (20 GeV/c) as measured by other experiments [16,17]. Therefore, to present a fair comparison, all of the predictions are adjusted to include this contribution using the FONLL calculation [38], shown as a green band in Fig. 14 (a).   The bars and boxes indicate the statistical and total systematic uncertainty, respectively. A common luminosity uncertainties of 11.2% and 8.1% for the µ + µ − and e + e − decay channels are not included.

V. COMBINED RESULTS
The scaling behavior of particle production with x T = 2p T / √ s is characteristic for production through fragmentation due to hard scatterings. The x T scaling (E d 3 σ dp 3 = g(x T )/s n/2 ) has been tested for pions, protons, and the J/ψ for various collision energies [39], where n is a free parameter which can be interpreted as the number of active partons involved in hadron production. Figure 15 shows the x T dependence of protons, pions, and J/ψ. The J/ψ measured in 510 and 500 GeV proton+proton collisions has been fit to extract the parameter n, with n = 5.6 ± 0.1. This value is consistent with n = 5.6 ± 0.2 as found in a previous STAR measurement [40], as well as other previous measurements [32,[41][42][43][44][45][46][47][48][49] at high-p T , and this value is close to the CO and CEM predictions, which are n ∼ 6 [50,51] and smaller than that from NNLO ⋆ CSM prediction which is n ∼ 8 [52]. The broken scaling at low-p T is due to the onset of soft processes [40].

VI. CONCLUSIONS
Differential cross sections for the J/ψ meson in pro-ton+proton collisions at √ s = 510 and 500 GeV at RHIC are measured using the µ + µ − and e + e − decay channels. The results cover a wide p J/ψ T range from 0 to 20 GeV/c within |y J/ψ | < 0.4 and 1.0 for the µ + µ − and e + e − channel, respectively. Two different measurements of the inclusive J/ψ production cross section have been presented. The first is a fiducial cross section measurement utilizing only a restricted phase space as defined by detector acceptance, and is independent of the assumptions regarding polarization which results in a large systematic uncertainty at low p T . This allows direct comparisons between measurements and theoretical calculations in the future, for more discriminating tests of the models. The second is a full cross section measurement, accessing the full J/ψ decay phase space, depending highly on assumptions regarding polarization. The integrated fiducial and full production cross sections measured for inclusive J/ψ within 0 < p J/ψ T < 9 GeV/c are 10.3 ± 0.9 (stat.) ± 1.6 (sys.) ± 1.1 (lumi.) nb and 67 ± 6 (stat.) ± 10 (sys.) +200 −18 (pol.) ± 7 (lumi.) nb, respectively, via the µ + µ − channel. For 4 < p J/ψ T < 20 GeV/c, they are 2.90 ± 0.08 (stat.) ± 0.22 (sys.) ± 0.24 (lumi.) nb and 10.7 ± 0.5 (stat.) ± 0.8 (sys.) +5.9 −2.2 (pol.) ± 0.9 (lumi.) nb, respectively, via the e + e − channel. The calculations from CGC+NRQCD, NLO NRQCD and ICEM [5,8], which cover low, high, and both p T regions respectively, give a reasonable description for the data within the polarization envelope. The x T dependence for J/ψ is also presented for 510 and 500 GeV proton+proton collisions and the result is consistent with measurements at other collision energies from other collaborations. The ratio of ψ(2S) to J/ψ for p T from 4−12 GeV/c is measured to be 0.038 ± 0.010 (stat.) ± 0.003 (sys.). It is consistent with results from other experiments and there is no obvious collision energy dependence. Since the J/ψ production mechanism is not yet fully understood, it is important to continue confronting the models that incorporate the most current understanding with new data. The results presented in this paper, together with cross section measurements at other energies, and measurements of the polarization, contribute to the goal of better understanding the production of heavy quarkonium.