Polarization and cross section of midrapidity J/$\psi$ production in proton-proton collisions at $\sqrt{s}=510$ GeV

The PHENIX experiment has measured the spin alignment for inclusive $J/\psi\rightarrow e^{+}e^{-}$ decays in $p$+$p$ collisions at $\sqrt{s}=510$ GeV at midrapidity. The angular distributions have been measured in three different polarization frames, and the three decay angular coefficients have been extracted in a full two-dimensional analysis. Previously, PHENIX saw large longitudinal net polarization at forward rapidity at the same collision energy. This analysis at midrapidity, complementary to the previous PHENIX results, sees no sizable polarization in the measured transverse momentum range of $0.0<p_T<10.0$ GeV/$c$. The results are consistent with a previous one-dimensional analysis at midrapidity at $\sqrt{s}=200$ GeV. The transverse-momentum-dependent cross section for midrapidity $J/\psi$ production has additionally been measured, and after comparison to world data we find a simple logarithmic dependence of the cross section on $\sqrt{s}$.

The PHENIX experiment has measured the spin alignment for inclusive J/ψ → e + e − decays in proton-proton collisions at √ s = 510 GeV at midrapidity. The angular distributions have been measured in three different polarization frames, and the three decay angular coefficients have been extracted in a full two-dimensional analysis. Previously, PHENIX saw large longitudinal net polarization at forward rapidity at the same collision energy. This analysis at midrapidity, complementary to the previous PHENIX results, sees no sizable polarization in the measured transverse momentum range of 0.0 < pT < 10.0 GeV/c. The results are consistent with a previous one-dimensional analysis at midrapidity at √ s = 200 GeV. The transverse-momentum-dependent cross section for midrapidity J/ψ production has additionally been measured, and after comparison to world data we find a simple logarithmic dependence of the cross section on √ s.

I. INTRODUCTION
Measurements of heavy quark bound states in hadronic collisions provide unique tools for testing Quantum Chromodynamics (QCD). Charmonium, the bound state of a charm and anticharm quark, is produced predominantly via gluon fusion at Relativistic Heavy Ion Collider (RHIC) kinematics. The J/ψ is a colorless neutral meson with spin 1 and decays with considerable branching ratio into lepton pairs (≈ 6% each for dielectrons and dimuons). Various theoretical models have been developed to describe the J/ψ production cross section and polarization as a function of transverse momentum, but none can describe both simultaneously; see Refs. [1] and [2] for recent reviews of theoretical developments and phenomenological work.
One theoretical approach to J/ψ production is within the framework of nonrelativistic QCD (NRQCD), which is an effective theory that describes the dynamics of heavy quark bound states at nonrelativistic scales (ν = v/c 1) [3][4][5][6]. The large heavy-quark mass scale relative to the hadronization scale, m Q Λ QCD , factorizes the J/ψ production process into the quark-antiquark (QQ) pair production at short-distance and subsequent formation of the heavy quark meson at long-distance. In the former regime, process-dependent short-distance coefficients are calculated perturbatively, and in the latter nonperturbative regime, the behavior is encoded in longdistance matrix elements (LDMEs). The QQ intermediate states are allowed to have quantum numbers different from those of the final-state meson. The leading-order relativistic corrections, for instance, put the QQ either in the color-singlet (CS) state 3 S    J (J=0,1,2). The relative importance of states with different color and angular momentum quantum numbers is estimated by νscaling rules [4]. Including CO states is found to be cru- * Deceased † PHENIX Spokesperson: akiba@rcf.rhic.bnl.gov cial as their leading-power contributions ≈ O( 1 p T ) 4 show up at leading order (LO), in contrast to the CS terms that only appear at next-to-next-to-leading order (NNLO). The LDMEs can only be determined from experimental data by fitting them via a global analysis. Several groups who performed global analyses [7][8][9][10] successfully described world data for transverse momentum (p T ) spectra, while consistent predictions of the spin alignment measurements for quarkonia remain a challenge.
The spin alignment of a positively charged lepton from a J/ψ decay, commonly known as "polarization," has been measured at the Tevatron [11], RHIC [12,13], and the Large Hadron Collider [14][15][16]. Measuring spin alignment provides additional tests for the theory and understanding dominant quarkonium production mechanisms in different kinematic regimes. The J/ψ polarization is measured by fitting the angular distribution of a positively charged lepton, shown in Eq. (1), to data and extracting decay angular coefficients. dN dΩ ≈ 1 + λ θ cos 2 θ + λ θφ sin 2 θ cos 2φ + λ φ sin 2θ cos φ, (1) where the coefficients λ θ , λ θφ , and λ φ are determined most commonly in the helicity (HX) frame [17], Collins-Soper (C-S) frame [18] and Gottfried-Jackson (G-J) frame [19] defined in the J/ψ production plane. Invariant variables are constructed using SO(2) symmetry in choosing the z-axis in the production plane. Physical interpretation is straightforward only with these invariant quantities, making direct comparison possible between experimental results in different kinematic regimes. Equation (2) shows the two frame-invariant variables defined in Refs. [20,21] to characterize the decay angular distribution.
Recently, general methods have been developed for finding all independent invariants under rotations for particles with various spin quantum numbers [22,23]. There are also new theoretical developments for less inclusive observables that include looking at the polarization of quarkonia produced in jets [24][25][26][27].
Previously a PHENIX J/ψ measurement at √ s = 510 GeV in p+p at forward rapidity showed largely longitudinal net polarization (negative angular coefficients) [13], while a prior midrapidity PHENIX measurement at √ s = 200 GeV was consistent with no strong polarization [12]. The present analysis for midrapidity J/ψ production in √ s = 510 GeV collisions is complementary to both previous measurements.

II. EXPERIMENTAL SETUP
In 2013, the PHENIX experiment at RHIC collected data from longitudinally polarized p+p collisions at √ s = 510 GeV. An integrated luminosity of 136 pb −1 was used for J/ψ polarization measurements at midrapidity. The PHENIX detector is described in detail in Ref. [28]. The PHENIX central arms comprise the east and west arms and cover a pseudorapidity range |η| < 0.35 and an azimuthal coverage of ∆φ = π 2 for each arm. The PHENIX detector elements used in this analysis include the drift chamber (DC), pad chambers (PC), ring-imagingČerenkov (RICH) detector and electromagnetic calorimeter (EMCal). Charged particle tracks are reconstructed with the DC and PC tracking system. These detectors also provide the momentum information of the tracks. The data sample used for this analysis is the sum of events triggered with three different energy deposit thresholds applied in the EMCal. Triggers used in this study are described in more detail in Section III and in Ref. [29]. The RICH was used for electron identification. Two sets of 64 quartz-crystal radiators attached to photomultipliers at z positions of ± 144 cm and pseudorapidity 3.1 < |η| < 3.9 were used to trigger hard collision events and to evaluate the collision vertex position. These beam-beam-counters (BBCs) and zerodegree calorimeters (ZDCs) were used together to evaluate and compare the luminosities seen by PHENIX. In addition, the silicon vertex detector (VTX) was placed in the west arm around the beam pipe at nominal radii of 2.6, 5.1, 11.8, 16.7 cm with an acceptance of |η| < 1 and ∆φ = 0.8π. The total material budget expressed as a percentage of a radiation length is 13.42%. As the VTX was not in operation in 2013, this created a large source of electron background from conversions of direct and decay photons.

A. Event Selection
Events were triggered by particles leaving energy of at least 2.2, 3.7 or 4.7 GeV in the EMCal for the various EMCal-RICH trigger (ERT) thresholds. Different scale-down factors were used for each threshold, that is, different fractions of triggered events were written to tape. An OR of all these triggers was used for polarization measurements and only the lowest threshold trigger was required for cross section measurements. The J/ψ candidates are triggered by either or both members of the electron-positron pair. The transverse momentum of electrons is determined by the bending of the track in the magnetic field before the DC. For electron identification, at least one photomultiplier tube in the RICH, associated with the track, is required to have fired.
Due to the partial azimuthal coverage on each side of the PHENIX Central Arm detector, different arm combinations in the measurements of decay pair products give access to different pair p T regions. For instance, high-p T enhancement is achieved in the sample of J/ψ mesons whose decay dielectrons are detected in the same arm due to their small opening angle in the lab frame. The dielectrons detected in the east arm offer a relatively clean high-p T -enhanced J/ψ sample. The pairs within the west arm were not used in this analysis due to excessive background consisting of γ to e + e − conversion caused by presence of the VTX that was not operational. The pairs with one lepton in each arm were used to obtain J/ψ mesons with lower p T .
The number of J/ψ events is determined by fitting the invariant mass of lepton pairs in the data and counting within the mass range between 2.8 and 3.3 GeV/c 2 . Combinatorial background attributed to decay electrons from hadrons is eliminated statistically before performing the fit using the like-sign method [30] established in PHENIX. The sum of the signal and residual background components, described by a Crystal Ball function [31] and an exponential function, respectively, is used to further subtract correlated electron background coming from open heavy flavor and Drell-Yan decay processes. Invariant mass distributions for dielectrons measured in the same arm and opposite arms are shown in Fig. 1 along with fit results.

B. Efficiency and acceptance corrections
The MC pseudodata were generated to estimate the effects of detector acceptance and efficiency. Events containing a single J/ψ per event were generated with zero polarization, a flat p T distribution, and a flat z-vertex distribution along the beam direction, using Monte Carlo techniques. The J/ψ's subsequently were forced to decay into dielectrons within the detector acceptance. The yields were then reweighted, as described in more detail below, to emulate effects that are present in the data such as smeared vertex distribution along the beam axis, the shape of the J/ψ p T distribution and the ERT trigger efficiencies.
Polarization measurements are particularly susceptible to inconsistency between the simulation and data as the acceptance corrections accounting for the zero polarization baseline are solely dependent on the simulation. A comparison of the polar and azimuthal angular distributions between the data and MC simulation, after the effects of dead areas, shows good agreement, as seen in Fig. 2.
The shape of dN (J/ψ) dp T can affect the polarization extraction in a nontrivial manner due to the limited azimuthal acceptance of the PHENIX Central Arms. Opposite-sign pairs bending toward each other into the detector populate different phase space than the ones bending away from each other. For this reason, the p T differential J/ψ cross section has been measured in data. The measured cross section was fit with a Kaplan function, defined in Eq. 3, and the fit results were then used to reweight the MC pseudodata for the acceptance correction.
As the data were sampled from an OR of three triggers with different energy thresholds, a method was developed to emulate the trigger efficiency effect on the p T shape of the data. The efficiency of triggering on J/ψ events for each type of ERT was determined from the minimum-bias dataset by checking the trigger condition. The acceptance-corrected pseudodata were then processed by randomly sampling the distributions of ERT-measured efficiencies and prescale factors. The detector area masked due to being ineffective or suffering from excessive electron background was taken into account in this method.
The shape of the raw J/ψ yield as a function of p T in data before corrections is in excellent agreement with the efficiency corrected normalized yield in simulation for each trigger type as well as for all possible trigger combinations (see Fig. 3).

C. Cross Sections
For cross section measurements, the raw J/ψ yield is corrected for trigger efficiencies as well as track reconstruction efficiencies and acceptance. The EMCal trigger efficiencies were measured using data and the electron reconstruction efficiencies and acceptance corrections were obtained from MC simulation as described in the previous section. An integrated luminosity of 136 pb −1 was sampled in the analysis and the global normalization of the cross section was determined by the p+p cross section sampled by the BBC trigger, which is found to be 32.5 ± 3.0 (stat) ± 1.2 (syst) mb, based on van-der-Meer scan results. The effects of multiple collisions per beam crossing due to high luminosity at √ s = 510 GeV, which were estimated to be of order of 20%, were also taken into account in counting J/ψ yields.

D. Angular coefficients
To account for acceptance and efficiency effects, the raw angular distribution in data is divided by the reweighted single J/ψ MC pseudodata generated under an assumption of no polarization as described in Section III B. This procedure closely follows methods de-scribed in Ref. [13] and additionally adopts the trigger emulator described in Section III B. The cos θ − φ distributions of the decay positron in three different polarization frames in MC pseudodata are shown in Fig. 4.
Angular coefficients were extracted with three different methods; the χ 2 method, the maximum loglikelihood (MLL) method and the parametric bootstrap method [32]. The χ 2 fitting results are displayed in contours at 68.3%, 95.4% and 99.7% confidence level (CL) in Fig. 5. One can see that some of the coefficients are correlated in one frame and uncorrelated in other frames. The helicity frame is known to be orthogonal to the C-S frame and similar to the G-J frame. This is also seen in our results, i.e. the correlation pattern between λ θ and λ φ , represented by the orientation of elliptic fit contours, in the helicity frame appears to be better aligned with the ones in the G-J frame. The MLL method treats low statistics bins more properly than the χ 2 method, while it is prone to underestimating uncertainties as the signalto-background ratio becomes smaller. In fitting methods, the uncertainties on the invariants that are a function of primary angular coefficients are propagated from the ones on coefficients, taking into account the correlation matrix. The parametric bootstrap method properly estimates uncertainties on the invariant variables. In this method, a Gaussian noise term is added to each measurement point at each sampling trial with a standard deviation of the corresponding measured uncertainty. After a large number of sampling procedures, the uncertainties are defined at a 68% CL.
The results are consistent among different methods within 0.5σ in all frames and the central measured values and statistical uncertainties from the best-fit results were taken as the final results.
At lower p T (<3.0 GeV/c), the decay electron and A silicon vertex detector (VTX) was installed in the PHENIX west arm, but was not operational during the data taking. This resulted in a high level of irreducible electron background originating from photon conversion in the VTX material. These random combinations of background electrons were subtracted using an event mixing technique, but resulted in higher statistical and additional systematic uncertainty for the low p T measurement. The PHENIX acceptance in cos θ − φ space is very different for decay electrons and positrons detected in opposite PHENIX Central Arms compared to the case when both the electron and positron are detected in the same arm. As a result, the acceptance for the C-S frame turned out to be very limited, and a polarization measurement in this frame was not done at low p T . The polarization coefficients at low p T were determined using only the χ 2 fitting method.

E. Systematic uncertainties
Sources of systematic uncertainties on the polarization and cross section measurement include uncertainties on acceptance, tracking and electron identification efficiency, trigger efficiencies, uncertainties on the global normalization derived from the cross section seen by the BBC, and the trigger rate dependence of J/ψ yield counting. The two latter sources do not affect the polarization measurements at all, but they are the main sources of systematic uncertainty for J/ψ cross section determination.
During RHIC run 2013, the proton beam luminosity was very high, which resulted in a nonnegligible probability of multiple collisions per beam crossing at high trigger rates. The correction to the J/ψ cross section due to multiple collisions per crossing was calculated by dividing the whole data set into five run groups with different trigger rates, calculating the cross section for each group, and extrapolating to zero trigger rate. The correction was estimated to be 0.80 ± 0.20. The large uncertainty of this estimate is due to extrapolating from very high trigger rates to zero and lack of statistics. This is the largest systematic uncertainty of the J/ψ cross section.
The primary source of systematic uncertainties on the polarization measurements in the central-arm detectors is the shape of the J/ψ p T spectra. The p T shape is af-fected by statistical uncertainties on the raw measured yields as well as uncertainties on the ERT efficiencies. The uncertainties on the ERT efficiencies were estimated using a sampling method. In each sampling trial, parameters were extracted from fitting the efficiency curve for each trigger type, and the analysis was repeated with those parameters to estimate the combined ERT efficiencies. The quadratic sum of the statistical and systematic uncertainties was assigned as a total uncertainty at each measurement point. The systematic uncertainties on angular coefficients that are attributed to p T shape were estimated using the parametric bootstrap method with a Gaussian noise term corresponding to the aforementioned uncertainties. Uncertainties estimated in this method depend on the p T of J/ψ, the type of λ coefficient and the polarization frame. Resulting uncertainties from this source are larger (as high as ≈0.5 on λ θ ) at p T < 3.0 GeV/c and <0.1 in most cases at p T ≥ 3.0 GeV/c. In addition, the differences between fitting methods, estimated to be less than 0.5 times the statistical uncertainties in all angular coefficients, were added as systematic uncertainties at high p T .
At low p T (below ≈3.0 GeV/c), the PHENIX acceptance in cos θ − φ space is very different from that at higher p T . Relative contributions from different sources of systematic uncertainties also change, while the dependence on exact p T shape remains a dominant factor. An additional systematic uncertainty at low p T is caused by the much larger combinatorial background, which had to be subtracted using a mixed-event method. This uncertainty was estimated by changing the mixed-event background normalization and was found smaller than the two uncertainties mentioned above. Systematic uncertainties of this measurement are shown in Table I.

IV. RESULTS
The measured p T -differential J/ψ cross section is shown in Fig. 6 with a fit to a Kaplan function. Due to the higher collision energy the p T spectrum is harder than previously published PHENIX results for √ s = 200 GeV [12], where the fit parameter b that determines the hardness of the spectrum for a given n was estimated to be smaller at 3.41 ± 0.21 and the parameter n was comparable at 4.6 ± 0.4. Figure 6 also shows a comparison of the measured J/ψ differential cross section with a theory prediction based on full NRQCD at NLO with leading relativistic corrections that includes CS and CO states, provided by Butenschön et al. [33]. Within its uncertainties, the theory calculation is in agreement with the experimental results within its valid range of p T > ∼ 2 GeV/c, justifying use of the theory model for predictions of polarization measurements in this kinematic range.
The p T -integrated cross section times branching ratio is shown in Fig. 7 along with the previous PHENIX results at √ s = 200 GeV [34] and the world results from the Large Hadron Collider [35,36] and Tevatron [37]. A sim-ple logarithmic dependence on the collision energy is seen for J/ψ production at midrapidity, making estimates of J/ψ yield at any √ s easy, and inviting the theory community to model the trend. The p T -integrated cross section times branching ratio was found to be 97.6 ± 3.6 (stat) ± 5.1 (syst) ± 9.8 (global) ± 19.5 (multiple collision) nb.
FIG. 6. The measured differential J/ψ cross section times branching ratio as a function of transverse momentum. Fit parameters are A = 37.6 ± 2.2 nb/(GeV/c), b = 4.33 ± 0.28 GeV/c, and n = 4.61 ± 0.32. Open rectangles show systematic point-by-point uncertainties. Global normalization uncertainty of 10.1% is not shown. Green curve shows theory prediction based on full NRQCD at NLO with leading relativistic corrections that include CS and CO states, provided by Butenschön et al. [33]. Light green band indicates theory uncertainty.
The final results of the three primary angular coefficients are shown in Fig. 8, Fig. 9, and Fig. 10. Uncertainties on measurements in different frames are correlated. Due to limited detector acceptance in η at midrapidity, λ θ (Fig. 8) is poorly constrained compared to λ φ (Fig. 9). Theory predictions based on full NRQCD at NLO with leading relativistic corrections that includes CS and CO states, provided by Butenschön et al. [38], are overlaid with the measurements. The uncertainty bands on the NRQCD predictions account for the scale uncertainties and uncertainties on the LDMEs. The LDMEs were obtained in a global analysis of unpolarized data that excludes measurements from hadroproduction with p T <3.0 GeV/c. To improve consistency with the data, the feed-down from heavier charmonium states was subtracted in the theory prediction. The fraction of J/ψ events from b-flavored hadron decays is negligible at RHIC. In PHENIX, the unpolarized yield measurements [12] are well described down to 1 GeV/c, justifying     Measuring these invariant variables not only provides a direct test for the underlying production mechanisms, but also robust tools to address systematic uncertainties by comparing results in different frames. Consistency between different frames shown forλ in Fig. 11 and for F in Fig. 12 indicates that systematic effects are under good control. Frame-invariant variables are interpreted as the total net polarization summed over all production mechanisms. Different production mechanisms can lead to natural polarization in different frames. In the case ofλ, transverse polarization corresponds to the value of 1 and longitudinal polarization to the value of -1. The zero value seen inλ is interpreted as no net polarization. The other variable F also carries similar meaning. Unlikẽ λ, F is mathematically bounded between 0 and 1. The zero polarization corresponds to 1 3 and transverse and longitudinal polarization to 2 3 and 0, respectively. Two scenarios were considered for these frameinvariant results, NRQCD and the Color Singlet Model (CSM), which is equivalent to the v → 0 limit of NRQCD. The uncertainties on CSM predictions include the scale uncertainty. The CSM prediction is qualitatively opposite to the NRQCD predictions [38,39]. NRQCD predicts transverse polarization while the CSM predicts longitudinal polarization, as is shown in Fig. 11 and Fig. 12. Neither strong transverse nor longitudinal polarization is measured in the data in any frame. Results are also compared to previous PHENIX measurements at a different beam energy [12] and rapidity [13]. Figure 13 compares λ θ in the helicity frame for different collision energies. At midrapidity, the de-FIG. 14. Angular coefficient λ θ of J/ψ production at √ s = 510 GeV shown as a function of rapidity in three polarization frames. The points for different frames are shifted for clarity. Forward rapidity points are from Ref. [13].
FIG. 15. Rapidity dependence ofλ of J/ψ production at √ s = 510 GeV in three polarization frames. The points for different frames are shifted for clarity. Forward rapidity points are from Ref. [13]. Frame pT λ θ λ φ λ θφλ F Helicity 0.0< pT <3.0 0.01 ization both at √ s = 200 GeV (1-dimensional analysis) and 510 GeV (2-dimensional analysis). The rapidity dependence of λ θ andλ is shown in Fig. 14 and Fig. 15, respectively. While midrapidity data indicate no polarization, moderate polarization is seen at forward rapidity. At forward rapidityλ was measured to be largely negative, indicating longitudinal polarization. This is in stark contrast to this result that sees no preferred polarization direction, shown in Fig. 14 and Fig. 15. No strong polarization was seen in other experiments at higher p T and higher beam energies in general, and the discrepancy between measurements and theory predictions is still being studied. Results of polarization measurements are summarized in Table I.

V. SUMMARY
The PHENIX experiment measured the J/ψ polarization at midrapidity in √ s = 510 GeV p+p collisions by reconstructing the hadronized charmonium state in the dielectron decay channel. The midrapidity cross section at √ s = 510 GeV in the same channel has been newly measured and is consistent with NRQCD calculations above p T > ∼ 2 GeV/c. The results show the expected hardening of the J/ψ p T spectrum as compared to the measurement at √ s = 200 GeV. At both low p T and high p T , the net polarization has been observed to be consistent with zero within uncertainties. This is in contrast to the measurements made at forward rapidity. The new results do not rule out either the CSM or the NRQCD J/ψ production models. The new measurements from the 2-dimensional analysis show consistency in λ θ with the results from a previous 1-dimensional midrapidity analysis at √ s = 200 GeV.