$J/\psi$ and $\psi(2S)$ production at forward rapidity in $p$+$p$ collisions at $\sqrt{s}=510$ GeV

The PHENIX experiment at the Relativistic Heavy Ion Collider has measured the differential cross section, mean transverse momentum, mean transverse momentum squared of inclusive $J/\psi$ and cross-section ratio of $\psi(2S)$ to $J/\psi$ at forward rapidity in \pp collisions at \sqrts = 510 GeV via the dimuon decay channel. Comparison is made to inclusive $J/\psi$ cross sections measured at \sqrts = 200 GeV and 2.76--13 TeV. The result is also compared to leading-order nonrelativistic QCD calculations coupled to a color-glass-condensate description of the low-$x$ gluons in the proton at low transverse momentum ($p_T$) and to next-to-leading order nonrelativistic QCD calculations for the rest of the $p_T$ range. These calculations overestimate the data at low $p_T$. While consistent with the data within uncertainties above $\approx3$ GeV/$c$, the calculations are systematically below the data. The total cross section times the branching ratio is BR $d\sigma^{J/\psi}_{pp}/dy (1.2<|y|<2.2, 0

The PHENIX experiment at the Relativistic Heavy Ion Collider has measured the differential cross section, mean transverse momentum, mean transverse momentum squared of inclusive J/ψ and cross-section ratio of ψ(2S) to J/ψ at forward rapidity in p+p collisions at √ s = 510 GeV via the dimuon decay channel. Comparison is made to inclusive J/ψ cross sections measured at √ s = 200 GeV and 2.76-13 TeV. The result is also compared to leading-order nonrelativistic QCD calculations coupled to a color-glass-condensate description of the low-x gluons in the proton at low transverse momentum (pT ) and to next-to-leading order nonrelativistic QCD calculations for the rest of the pT range. These calculations overestimate the data at low pT . While consistent with the data within uncertainties above ≈ 3 GeV/c, the calculations are systematically below the data. The total cross section times the branching ratio is BR dσ J/ψ pp /dy(1.2 < |y| < 2.2, 0 < pT < 10 GeV/c) = 54.3 ± 0.5 (stat) ± 5.5 (syst) nb.

I. INTRODUCTION
Charmonium states such as J/ψ and ψ(2S) mesons are bound states of a charm and anti-charm quark (cc). At the Relativistic Heavy Ion Collider (RHIC) energies, they are produced mostly from hard scattering of two gluons into a cc pair followed by the evolution of this pair through a hadronization process to form a physical charmonium. Despite several decades of extensive studies [1][2][3][4][5][6][7][8][9] since the discovery of J/ψ, we still have very limited knowledge about the J/ψ production mechanism and hadronization. Therefore, carrying out as many charmonium measurements as possible in p+p collisions over a wide range of transverse momentum (p T ) and of rapidity (y) at different energies is essential to understanding production mechanisms. These measurements over a wide range of p T (down to zero p T ) and rapidity allow calculating quantities, such as the mean transverse momentum p T , the mean transverse momentum squared p 2 t , and the p T -integrated cross section dσ/dy. The collision energy dependence of these quantities can put stringent constraints on the different theoretical approaches that are used to describe the hadronic production of J/ψ. These approaches include the color-evaporation model (CEM) [10,11], the color-singlet model (CSM) [12] and the nonrelativistic quantum chromodynamics formalism (NRQCD) [13]. In this work, we compare the data to NRQCD, an effective field theory derived from QCD and valid for heavy-quark pairs with low relative velocity, where a J/ψ can be formed from cc pair produced in a color-singlet or a color-octet state.
In this paper, we present the inclusive J/ψ production cross section and the ratio of ψ(2S) to J/ψ production cross sections at forward rapidity (1.2 < |y| < 2.2) measured in p+p collisions at center of mass energy √ s = * Deceased † PHENIX Spokesperson: akiba@rcf.rhic.bnl.gov 510 GeV. These mesons are measured in the dimuon decay channel. The J/ψ inclusive differential cross sections are obtained as a function of p T and y over a wide range of p T . The J/ψ and ψ(2S) results at √ s = 510 GeV are the first measurements at this rapidity. Comparisons to similar PHENIX measurements performed at √ s = 200 GeV [2] and Large Hadron Collider (LHC) measurements at √ s = 2.76, 5.02, 7, 8 and 13 TeV [3][4][5][6] allow studying the variations of p T , p 2 t and dσ/dy as a function of √ s. The results are also compared to next-to-leading order (NLO) NRQCD calculations [8].
The paper is organized as follows: the PHENIX apparatus is described in Sec. II, the data samples used for this analysis and the analysis procedure are discussed in Sec. III, while the results are presented and compared to measurements at different √ s as well as to models in Sec. IV.

II. EXPERIMENTAL SETUP
A complete description of the PHENIX detector can be found in Ref. [14]. Only the detector systems relevant to this measurement are briefly described here.
The PHENIX muon spectrometers, see Fig. 1, cover the full aziumth and the north (south) arm cover forward (backward) rapidity, 1.2 < y < 2.2 (−2.2 < y < −1.2). Each muon spectrometer comprises a hadronic absorber, a magnet, a muon tracker (MuTr), and a muon identifier (MuID). The absorbers comprise layers of copper, iron and stainless steel and have about 7.2 interaction lengths. Following the absorber in each muon arm is the MuTr, which comprises three stations of cathode strip chambers in a radial magnetic field with an integrated bending power of 0.8 T·m. The MuID comprises five alternating steel absorbers and Iarocci tubes. The composite momentum resolution, δp/p, of particles in the analyzed momentum range is about 5%, independent of momentum and dominated by multiple scattering. Muon candidates are identified by reconstructed tracks in the MuTr matched to MuID tracks that penetrate through to the last MuID plane.
Since 2012 the PHENIX detector had a new forward vertex detector (FVTX) [15], which comprises four planes of silicon strip detectors, finely segmented in radius and coarsely segmented in azimuth. For the subset of muon candidate tracks passing several of these detector planes, this additional information was used to improve mass resolution by a factor of 1.5 for studying ψ(2S).
Another detector system relevant to this analysis is the beam-beam counter (BBC), comprising two arrays of 64Čerenkov counters, located on both sides of the interaction point and covering the pseudorapidity range 3.1 < |η| < 3.9. The BBC system was used to measure the p+p collision vertex position along the beam axis (z vtx ), with 2 cm resolution, and initial collision time. It was also used to measure the beam luminosity and form a minimum bias (MB) trigger.

III. DATA ANALYSIS
The results presented here are based on the data sample collected by PHENIX during the 2013 p+p run at √ s = 510 GeV. The BBC counters provided the MB trigger, which required at least one hit in each of the BBCs. Events, in coincidence with the MB trigger, containing a muon pair within the acceptance of the spectrometer are selected by the level-1 dimuon trigger (MuIDLL1-2D) requiring that at least two tracks penetrate through the MuID to its last layer. The data sample, used in this analysis, corresponds to 3.02 × 10 12 MB events or to an integrated luminosity of 94.4 pb −1 .

A. Raw yield extraction
A set of quality cuts is applied to the data to select good p+p events and good muon candidates as well as to improve the signal-to-background ratio. Good p+p events are selected by requiring that the collision occurs in the fiducial interaction region |z vtx | < 30 cm as measured by the BBC. Each reconstructed muon candidate comprises a combination of reconstructed muon tracks in the MuTr and in the MuID. The MuTr track is required to have more than 9 hits out of the maximum possible of 16 while the MuID track is required to have more than 6 hits out of the maximum possible of 10. In addition, a cut on individual MuTr track χ 2 of 23 is applied. The MuTr track χ 2 is calculated from the difference between the measured hit positions of the track and the subsequent fit for each MuTr track. The MuTr tracks are then matched to the MuID tracks at the first MuID layer by applying cuts on maximum position and angle differences. Furthermore, there is a minimum allowed single muon momentum along the beam axis, p z , which is reconstructed and energy-loss corrected at the collision vertex, of 3.0 GeV/c corresponding to the momentum cut effectively imposed by the absorbers. Finally, a cut on the χ 2 of the fit of the two muon tracks to the common vertex of the two candidate tracks near the interaction point was applied.
The invariant mass distribution is formed by combining muon candidate tracks of opposite charges (unlikesign). In addition to the charmonium signal, the resulting unlike-sign dimuon spectrum includes correlated and uncorrelated pairs. In the J/ψ and ψ(2S) region, the correlated pairs arise from correlated semi-muonic decays of charmed hadrons, beauty and the Drell-Yan process as well as light hadron decays. The uncorrelated pairs are mainly coming from the decays of light hadrons (π ± , K ± and K 0 ) which decay before or after passing through the absorber, and form the combinatorial background.
The combinatorial background is estimated using two methods: The first one derives the combinatorial background from the mass distribution of the same sign (like-sign) pairs of muon candidates within the same event. The second method derives the combinatorial background from the mass distribution of the unlike-sign pairs of muon candidates from different events (mixedevent) of z-vertex position within 2 cm. The normalization of the mass distribution of the combinatorial background from the like-sign dimuon distributions (N ++ and N −− ) is calculated as: N CB = 2 N ++ N −− . The mixedevent like-sign dimuon mass distribution is normalized to the same-event like-sign combinatorial background distribution in the invariant mass range 2.0−4.5 GeV/c 2 . This factor is then used to normalize the mixed-event unlikesign dimuon mass distribution. Figure 2 shows the unlike-sign dimuon spectrum together with the combinatorial background estimated by both methods. Both background distributions from the mixed-event and like-sign methods are consistent, however, the mixed-event background is more statistically stable, because we mix each event with the previous four events. Therefore, the mixed-event background was used to subtract the uncorrelated background from the unlikesign dimuon spectrum. After subtracting the uncorrelated background, the unlike-sign spectra including the correlated background are fitted by the following function, where p 0 − p 7 are free parameters and m µµ is the unlikesign dimuon mass. The J/ψ shape is better described with two Gaussian distributions, corresponding to the first two terms in Eq. 1, one for the J/ψ peak and a second one with larger width to account for the wider tails, which occurs due to limitations in MuTr resolution, as discussed in sec. II. The peak also includes contribution from ψ(2S), which is not resolved. An exponential is used to describe the continuum contributions from correlated backgrounds. Panels (a) and (b) of Fig. 2 show the raw spectra for selected p T and rapidity bins and panels (c) and (d) show the spectra after subtracting the combinatorial background fitted with the function described above for those selected bins.
To extract the ψ(2S) signal we improve the mass resolution of the muon tracking systems by utilizing the FVTX. The FVTX being located before the absorber allows measuring the dimuon opening angle before any multiple scattering occurs in the absorber [15]. Using this additional tracking information gives a more precise measurement of the dimuon opening angle and thereby a more precise measurement of the pair mass, as well as rejecting backgrounds from decay muons that emerge from the absorber. However, these additional requirements on the dimuon tracks that are necessary to separate the J/ψ and ψ(2S) peaks also reduce the statistics by a factor of 6 due to the geometric acceptance of FVTX, therefore, we study the dimuon mass spectra in each arm integrated over p T and rapidity within each arm. The dimuon mass spectrum extracted including the FVTX after subtracting the mixed-event background is shown in Fig. 3. Given the resolution enhancement, the sum of a Gaussian and a crystal-ball function [16,17], rather than a double Gaussian, was used for each of J/ψ and ψ(2S) peaks to fit the dimuon mass spectrum. The ψ(2S) peak is expected to be wider than the J/ψ peak, due to the fact that the higher mass and harder p T spectrum of the ψ(2S) state will produce higher momentum decay muons which have larger uncertainty in their reconstructed momentum in the spectrometer due to a smaller bend in the magnetic field. By selecting only poorly reconstructed tracks, we found a J/ψ width of ≈ 200 MeV/c 2 , therefore, the width of the second Gaussian in the fit to the entire sample of tracks is set to 200 MeV/c 2 . The ratio of widths of the ψ(2S) to J/ψ is set to 1.15, following expectations of the performance of the muon tracking system [18]. The difference between the centroids of the ψ(2S) and J/ψ peaks is set to the Particle Data Group value of 589 MeV/c 2 [19]. The relative normalization of the second Gaussian is fixed to be the same for both resonances, as are the parameters for the crystal-ball line shape.

B. Detector acceptance and reconstruction efficiency
The acceptance and reconstruction efficiency (Aε rec ) of the muon spectrometers, including the MuID trigger efficiency, is determined by running a pythia 1 [20] generated J/ψ signal through a geant4-based full detector simulation [21] of PHENIX. The simulation tuned the detector response to a set of characteristics (dead and hot channel maps, gains, noise, etc.) that described the performance of each detector subsystem. The simulated vertex distribution was tuned to match that of the 2013 data. The simulated events are reconstructed in the same manner as the data and the same cuts are applied as in the real data analysis.  the MuTr and MuID systems and different amount of absorber material.
In the case of ψ(2S), we are interested in the ratio of its differential cross section to that of J/ψ, therefore, we extract the ratio of Aε rec for ψ(2S) and J/ψ with addition of the FVTX information in analyzing the simulation to match that of the data analysis. A factor of 0.77 (0.69) is applied to the ψ(2S)/J/ψ ratio extracted from the fit to the invariant mass spectrum to account for differences in acceptance, efficiency, and dimuon trigger efficiencies between the north (south) arm of the muon spectrometer.

C. Differential cross section
The differential cross section is evaluated according to the following relation: where N ψ is the extracted J/ψ or ψ(2S) yield in y and p T bins with ∆y and ∆p T widths, respectively. BR is the branching ratio where BR J/ψ→µ + µ − = (5.93 ± 0.06) × 10 −2 and BR ψ(2S)→µ + µ − = (7.9±0.9)×10 −3 [19]. Aε rec is the product of the acceptance and reconstruction efficiency. N BBC MB = 3.02 × 10 12 is the number of MB events and ε BBC = 0.91±0.04 is the efficiency of the MB trigger for events containing a hard scattering [22]. σ BBC is the PHENIX BBC cross section, 32.5 ± 3.2 mb at √ s = 510 GeV, which is determined from the van der Meer scan technique [23].

D. Systematic uncertainties
All systematic uncertainties are evaluated as standard deviations and are summarized in Tables I and II. They are divided into three categories based upon the effect each source has on the measured results: Type-A: Point-to-point uncorrelated uncertainties allow the data points to move independently with respect to one another and are added in quadrature with statistical uncertainties; however, no systematic uncertainties of this type are associated with this measurement.
Type-B: Point-to-point correlated uncertainties which allow the data points to move coherently within the quoted range to some degree. These systematic uncertainties include a 4% uncertainty from MuID tube efficiency and an 8.2% (2.8%) from MuTr overall efficiency for the north (south) arm. A 3.9% signal extraction uncertainty is assigned to account for the yield variations when using different functions, i.e., second, third and fourth order polynomials, to fit the correlated background and ≈ 3% uncertainty is assigned to account for the ψ(2S) contribution. The systematic uncertainty associated with Aε rec includes the uncertainty on the input p T and rapidity distributions which is extracted by varying these distributions over the range of the statistical uncertainty of the data, yielding 4.4% (5.0%) for the north (south) arm. Additional 11.2% (8.8%) systematic effect for the north (south) arm was also considered to account for the azimuthal angle distribution difference between data and simulation.
To be consistent with the real data analysis, a trigger emulator was used to match the level-1 dimuon trigger for the data. The efficiency of the trigger emulator was studied by applying it to the data and comparing the resulting mass spectrum to the mass spectrum when applying the level-1 dimuon trigger which resulted in a 1.5% (2%) uncertainty for the north (south) arm. Type-B systematic uncertainties are added in quadrature and amount to 16.0% (12.4%) for the north (south) arm. They are shown as shaded bands on the associated data points.
Type-C: An overall normalization uncertainty of 10% was assigned for the BBC cross section and efficiency uncertainties [24] that allow the data points to move together by a common multiplicative factor. In the measurement of the ψ(2S) to J/ψ ratio, most of the mentioned systematic uncertainties cancel out. However, the fit that was used to extract the yields is more complex and additional systematic uncertainties arose from the constraints applied during the fitting process.
A systematic uncertainty from constraining the normalization factor is determined by varying the mass range over which the factor is calculated and a 3% systematic uncertainty is assigned for both arms. Systematic uncertainty of 3% (7%) was assigned to the north (south) arm on the fit range by varying the range around the nominal values, 2-5 GeV/c 2 . The effect of constraining the second Gaussian peak width to 200 MeV/c 2 was studied by varying the width between 175 and 225 MeV/c 2 , resulting in a systematic uncertainty of 12% (10%) in the north (south) arm.
The systematic uncertainty component on Aε rec that survived the ratio amounts to 2.7% (4.1%) in the north (south) arm. The systematic uncertainties associated with the ratio measurement are summarized in Table II.

IV. RESULTS
The inclusive J/ψ differential cross section as a function of p T is calculated independently for each muon arm, then the results are combined using the best-linearunbiased-estimate method [25]. Results obtained using the two muon spectrometers are consistent within statistical uncertainties. The combined inclusive J/ψ differential cross section is shown in Fig. 5 and listed in Table III. The gray shaded bands represent the weighted average of the quadratic sum of type-B systematic uncertainties of the north and south arms, ≈ 10.1%. The average is weighted based on the statistical uncertainties of each arm.
The data points are corrected to account for the finite width of the analyzed p T bins [26]. We compare the data to inclusive J/ψ data at 200 GeV [2] which show similar p T dependence. At low p T , the data are compared to prompt J/ψ leading-order (LO) NRQCD calculations [8,13] coupled to a Color Glass Condensate (CGC) description of the low-x gluons in the proton [9]. For the rest of p T range, the data are compared to prompt J/ψ NLO NRQCD calculations [8,13]. The LO-NRQCD+CGC calculations overestimate the data at low p T . The NLO-NRQCD calculations underestimate the data at high p T , while to some extent, are consistent with the data at intermediate p T , 3-5 GeV/c. It is important to stress that the nonprompt J/ψ contribution (from excited charmonium states and from B-meson decays) is not included in these calculations. This is expected to be a significant contribution at high p T ; therefore, the addition of the nonprompt J/ψ contribution could account for the difference between the data and calculations [27][28][29]. The p T coverage down to zero p T allows the extraction of the p T -integrated cross section, BR dσ J/ψ pp /dy(1.2 < |y| < 2.2, 0 < p T < 10 GeV/c) = 54.3 ± 0.5 (stat) ± 5.5 (syst) nb.
Inclusive J/ψ differential cross section as a function of rapidity is listed in Table IV and shown in Fig. 6, which also includes PHENIX inclusive J/ψ data at 200 GeV [2] and NLO-NRQCD calculations [8]. The 510 GeV data show a similar rapidity dependence pattern to that of the 200 GeV data. NLO-NRQCD calculations overestimate the data, and this is consistent with what was observed in the case of p T -dependent differential cross section (see Fig. 5) because the y-dependent differential cross section is dominated by the low-p T region where NRQCD calculation overestimates the data.
To quantify the feed-down contribution of excited charmonium states, the ratio of the cross section of ψ(2s) to J/ψ, multiplied by their respective branching ratio to dimuons, is measured (R = 2.84±0.45%) and shown in Fig. 7. This ratio is compared with other p+p and p+A systems at different collision energies [17,[30][31][32][33][34][35][36][37][38]. The results are consistent with world data within uncertainties with no significant dependence on collision energy.
To better understand the shape of the p T spectrum for J/ψ at forward rapidity and quantify its hardening at √ s = 510 GeV, we calculate the corresponding mean transverse momentum p T and mean transverse momentum squared p 2 T . This is done by fitting the inclusive J/ψ p T -dependent differential cross sections with the following function [2,6]:    [7,17,[30][31][32][33][34][35][36][37][38]. The associated uncertainties are the quadrature sum of the statistical and systematic uncertainties.   The first error is statistical, and the second is the systematic uncertainty from the maximum shape deviation permitted by the type-B correlated errors. Figure 8 shows p T as a function of √ s from this measurement compared with results from 200 GeV PHENIX data at the same rapidity range [2], and results from AL-ICE at different √ s values and in the rapidity range, 2.5 < y < 4.0 [42]. This result follows the increasing pattern observed between PHENIX results at 200 GeV and ALICE results at 2.76-13 TeV. Figure 9 shows p 2 T as a function of √ s from this measurement compared with several other measurements [1,2,6,39,40,42,43]. Similar to p T , p 2 T from this measurement follows the increasing pattern versus √ s established by several sets of data over a wide range of energies. Below √ s of 2 TeV, the trend is qualitatively consistent with a linear fit of p 2 T versus the log of the center of mass energy from Ref. [2]. However, above √ s of 2 TeV, the ALICE data indicate p 2 T grows at an increased rate which is interpreted by authors of Ref. [6] as due to the fact that ALICE data sets have different p T ranges. The bottom cross section also increases with increasing √ s, changing the relative prompt and B-meson decay contributions to the inclusive J/ψ samples discussed here [27,44]. This may also contribute to the observed differences in the measured p 2 T .
The dσ J/ψ pp /dy measurement at √ s = 510 GeV offers an opportunity to test the center-of-mass energy dependence of the p T -integrated cross section. Moreover, it bridges the gap between RHIC data at 200 GeV and ALICE data starting at 2.76 TeV [3][4][5][6]. However, ALICE data are collected at mid (|y| < 0.9) and forward (2.5 < y < 4.0) rapidities and to have a proper comparison we interpolate the ALICE data to the PHENIX forward rapidity range, 1.2 < y < 2.2. This is done by fitting the pythia generated dσ/dy distribution at each energy to the data at the same energy with only the normalization as a free parameter. An example is shown in Fig. 10. We used several pythia [45] tunes including PHENIX default, tune-A, modified tune-A and atlas-csc [46]. After fitting each of these pythia tunes to the data, we extracted dσ/dy at 1.2 < y < 2.2, from these fits. The rms value of the extracted dσ/dy from the different fits is used in the comparison to RHIC data. The error on the rms value is the rms of the errors associated with the fit results.  Figure 11 shows that the data are well described by a power law, dσ

V. SUMMARY
We studied inclusive J/ψ production in p+p collisions at √ s = 510 GeV for 1.2 < |y| < 2.2 and 0 < p T < 10 GeV/c, through the dimuon decay channel. We measured inclusive J/ψ differential cross sections as a function of p T as well as a function of rapidity. The p T integrated differential cross section multiplied by J/ψ branching ratio to dimuons is BR dσ J/ψ pp /dy (1.2 < |y| < 2.2, 0 < p T < 10 GeV/c) = 54.3 ± 0.5 (stat) ± 5.5 (syst) nb. With these data measured over a wide p T range, we calculated p T , p 2 T and dσ/dy. The results were compared to similar quantities at different energies from RHIC and LHC to study their √ s dependence. These new measurements could put stringent constraints on J/ψ production models.
The inclusive J/ψ differential cross sections were compared to prompt J/ψ calculations. These calculations included LO-NRQCD+CGC at low p T and NLO-NRQCD for the rest of p T range. These model calculations overestimated the data at low p T and underestimated the data at high p T . The nonprompt J/ψ contribution was not included which could account for the underestimation at high p T where the nonprompt processes are significant.
In addition, we measured the ratio of the cross section of ψ(2s) to J/ψ, multiplied by their respective branching ratio to dimuons, R = 2.84 ± 0.45%. The result is consistent with world data within uncertainties with no dependence on collision energy.