Measurement of inclusive J/$\psi$ polarization in p+p collisions at $\sqrt{s}$ = 200 GeV by the STAR experiment

We report on new measurements of inclusive J/$\psi$ polarization at mid-rapidity in p+p collisions at $\sqrt{s}$ = 200 GeV by the STAR experiment at RHIC. The polarization parameters, $\lambda_\theta$, $\lambda_\phi$, and $\lambda_{\theta\phi}$, are measured as a function of transverse momentum ($p_T$) in both the Helicity and Collins-Soper (CS) reference frames within $p_T<10$ GeV/$C$. Except for $\lambda_\theta$ in the CS frame at the highest measured $p_T$, all three polarization parameters are consistent with 0 in both reference frames without any strong $p_T$ dependence. Several model calculations are compared with data, and the one using the Color Glass Condensate effective field theory coupled with non-relativistic QCD gives the best overall description of the experimental results, even though other models cannot be ruled out due to experimental uncertainties.


I. INTRODUCTION
The J/ψ meson, a bound state of a charm (c) and an anti-charm (c) quark, provides a natural testing ground for studying both the perturbative and non-perturbative aspects of the Quantum Chromodynamics (QCD). Due to their large masses, the production cross section of cc pairs can be calculated perturbatively. On the other hand, the formation of J/ψ mesons from cc pairs happens over long distances, and therefore is non-perturbative. The J/ψ mesons are also widely used in heavy-ion physics as an internal probe to study the properties of the quarkgluon plasma [1], which requires the measurement of the J/ψ production in vacuum as a reference. Despite decades of concentrated experimental and theoretical efforts, a complete picture of the J/ψ production mechanism in elementary collisions has yet to emerge.
Model calculations describing the J/ψ production utilize the factorization of the short-distance cc production and the long-distance hadronization process. Models differ mainly in the treatment of the non-perturbative formation of J/ψ. One of the early models is the Color Evaporation Model (CEM) [2,3], which is based on the principle of quark-hadron duality and satisfies all-order factorization. It assumes that every cc pair, with an invariant mass below twice the D-meson threshold, evolves into a J/ψ meson with a fixed probability (F J/ψ ) by randomly emitting or exchanging soft gluons with other color sources. The non-perturbative J/ψ formation is incorporated into the universal probability F J/ψ , which is independent of the kinematics and spin of the J/ψ meson. An Improved Color Evaporation Model (ICEM) has recently been proposed, in which the lower limit of the cc pair invariant mass is increased to be the charmonium mass and the transverse momentum (p T ) of the charmonium state is adjusted based on the ratio of its mass to the cc mass [4]. The ICEM calculation is in general agreement with the inclusive J/ψ cross section measured in p+p collisions at √ s = 200 GeV [4], in which the discrepancy seen above J/ψ p T ∼ 4 GeV/c is mainly due to the missing contribution of b-hadron decays in ICEM. By summing with the contribution of J/ψ from b-hadron decays obtained from the fixed-order plus next-to-leading logarithm (FONLL) calculation [5], the ICEM calculation agrees reasonably well with the inclusive J/ψ cross section measured in p+p collisions at √ s = 500 GeV [6] up to p J/ψ T = 20 GeV/c. A further extension based on ICEM at leading order (LO) is the calculation of J/ψ polarization utilizing the k T -factorization approach [7]. Compared to the measured J/ψ polarization at forward rapidity in p+p collisions at √ s = 7 TeV [8,9], the ICEM calculation shows significant discrepancies at low p T .
A more sophisticated way to describe the hadronization of heavy quarkonia is based on the effective quantum field theory of non-relativistic QCD (NRQCD) [10]. In addition to the usual expansion in the strong coupling constant (α s ), it also introduces an expansion in the relative velocity between the heavy quarks in the pair. Both the color-singlet and color-octet intermediate cc pairs are included in the NRQCD. The hadronization process is incorporated through the assumed universal Long Distance Matrix Elements (LDMEs), which weight the relative contributions of each intermediate state and are extracted from fitting experimental data. The NRQCD calculations at next-to-leading order (NLO) in α s have been done by three groups [11][12][13]. They obtained very different LDMEs depending on the low-p T cuts imposed on data points used and whether the polarization data are included. None of these calculations can give a simultaneous description of both the charmonium cross section, such as the η c yields measured in 7 TeV p+p collisions [14], and polarization such as those measured by the CDF Collaboration [15,16]. To remedy the issue of calculating the cc production cross section at low p T , where the collinear factorization formalism may not be applicable, an effort has been made to use the Color Glass Condensate (CGC) effective field theory [17]. Combined with the NRQCD, it describes well the J/ψ cross sections measured in p+p collisions at both RHIC and the LHC [18]. The CGC+NRQCD formalism has also been used to calculate the J/ψ polarization and the results agree well with the LHC measurements at forward rapidities [19]. Continued efforts from both experimental and theoretical sides are still needed to achieve the final goal of a complete understanding of J/ψ production.
While the J/ψ production cross section has been measured extensively in p+p collisions at √ s = 200 GeV at RHIC [20][21][22], its polarization, which is the topic of this paper, is less so [20,23]. The J/ψ polarization can be measured through the angular distribution of the positively charged daughter lepton [24]: (1 + λ θ cos 2 θ + λ φ sin 2 θ cos 2φ where λ θ , λ φ and λ θφ are the J/ψ polarization parameters. θ and φ are the polar and azimuthal angles of the positively charged daughter lepton in the J/ψ rest frame with respect to a chosen quantization axis. In the helicity (HX) frame [25], one usually uses the opposite of the direction of motion of the interaction point in the J/ψ rest frame as the quantization axis. In the Collins-Soper (CS) frame [26], one usually chooses the bisector of the angle formed by one beam direction and the opposite direction of the other beam in the J/ψ rest frame. J/ψ is considered fully transversely or longitudinally polarized when the polarization parameters take the values of (λ θ , λ φ , λ θφ ) = (1, 0, 0) or (-1, 0, 0). No polarization is referred to the case of (0, 0, 0). While the measured polarization values depend on the selection of the quantization axis, one can construct a frame invariant quantity to check the consistency of measurements in different frames [27]. It is defined as: Previous measurements of inclusive J/ψ polarization in 200 GeV p+p collisions [20,23] have only focused on λ θ in the HX frame within 0 < p T < 6 GeV/c. In this paper, we extend the scope by measuring all three polarization parameters in both HX and CS frames for p T < 10 GeV/c, as well as the frame invariant quantity λ inv . Measurements are carried out based on both the dimuon and dielectron decay channels covering different kinematic ranges. The inclusive J/ψ sample used in this paper includes directly produced J/ψ's and those from decays of excited charmonium states such as χ c and ψ(2S) (∼40% [28]) as well as b-hadrons (∼10-25% above p T of 5 GeV/c [21]).These results will provide more stringent tests of different model calculations, especially for the universality of model parameters, such as F J/ψ , LDMEs, that give models their predictive power. This paper is arranged as the following. An introduction to the Solenoidal Tracker At RHIC (STAR) is given in section II, followed by detailed descriptions of the analyses utilizing the electron and muon decay channels in sections III and IV, respectively. The J/ψ polarization results are presented in section V, and a summary is given in section VI.

II. STAR EXPERIMENT
The STAR experiment [29] at RHIC consists of a suite of mid-rapidity detectors with excellent tracking and particle identification (PID) capabilities. The Time Projection Chamber (TPC) [30] is the main tracking device to measure a particle's momentum and specific energy loss (dE/dx) for particle identification. It covers the pseudorapidity range of |η| < 1 over full azimuthal angle. A room temperature solenoidal magnet generates a uniform magnetic field of maximum value 0.5 T [31]. The Barrel Electromagnetic Calorimeter (BEMC) [32], covering full azimuthal angle over |η| < 1, is used to identify and trigger on high-p T electrons. In conjunction with the start time provided by the Vertex Position Detector (VPD), the Time-Of-Flight (TOF) detector [33] measures a particle's flight time to help further improve the electron purity. For the muon channel analysis, the Muon Telescope Detector (MTD) [34] is used for triggering on and identifying muons. It resides outside of the magnet which is used as an absorber, and covers about 45% in azimuth within |η| < 0.5. The forward-rapidity trigger detectors, the VPD at 4.24 < |η| < 5.1 [35] and Beam-Beam Counters (BBC) at 3.3 < |η| < 5.0 [36], are used to select collisions.

A. Dataset, event and track selections
The dataset was taken for p+p collisions at √ s = 200 GeV in 2012 using both the minimum-bias (MB) and high-tower (HT) triggers. The MB trigger selects non-single diffractive p+p collisions with a coincidence signal from the VPD on east and west sides, while the HT trigger selects events with energy depositions in the BEMC above given thresholds. About 300 million MB events, corresponding to an integrated luminosity of about 10 nb −1 , are analyzed to study the J/ψ polarization below p T of 2 GeV/c. Data collected by the HT0 (HT2) trigger with an energy threshold of E T > 2.6 (4.2) GeV correspond to an integrated luminosity of 1.36 (23.5) pb −1 . The HT0 trigger is used for the J/ψ measurement within 2 < p T < 4 GeV/c, while the HT2 trigger for 4 < p T < 14 GeV/c.
The vertex position along the beam direction can be reconstructed from TPC tracks (V TPC z ) or from the time difference of east and west VPD signals (V VPD z ). A cut of |V TPC z | < 50 cm is applied to ensure good TPC acceptance for all the events. An additional cut of |V TPC z − V VPD z | < 6 cm is applied to reduce the pile-up background from out-of-time collisions for MB events.
Charged tracks are required to have at least 20 TPC space points (out of a maximum of 45), a ratio of at least 0.52 between actually used and maximum possible number of TPC space points, at least 11 TPC space points for dE/dx calculation, and their distance of closest approach to the primary vertex (DCA) less than 1 cm. Electrons and positrons are identified using dE/dx in TPC, the velocity (β) calculated from the path length and time of flight between the collision vertex and TOF, and the ratio between the track momentum and energy deposition in the BEMC (pc/E). The normalized dE/dx is quantified as: where (dE/dx) measured is the measured energy loss in the TPC, (dE/dx) e theory the expected energy loss for an electron based on the Bichsel formalism [37], and σ(ln(dE/dx)) the resolution of the ln(dE/dx) measurement. The value of nσ e is required to be within (-1.9, 3). A cut of |1/β − 1| < 0.03 is applied for TOF-associated candidates, and 0.3 < pc/E < 1.5 for BEMC-associated candidates above 1 GeV/c. The electron and positron candidates are required to pass the nσ e cut, and either the β or pc/E cut. For HT-triggered events, at least one daughter of a J/ψ candidate must pass the pc/E requirement and have an energy deposition in the BEMC higher than the HT trigger threshold.

B. Analysis procedure
A maximum likelihood method is used to extract all three J/ψ polarization parameters simultaneously. The likelihood is defined as: where the sum is taken over the (cos θ, φ) bins, N J/ψ (cos θ, φ) is the raw number of J/ψ candidates in each (cos θ, φ) bin, and A × ε(cos θ, φ) the detector acceptance times J/ψ reconstruction efficiency in the same bin. F (cos θ, φ|λ θ , λ φ , λ θφ ) is the integral probability corresponding to cos θ and φ bin of positrons for given (λ θ , λ φ , λ θφ ) values, described by Eq. 1 normalized to 1.
A × ε(cos θ, φ) is evaluated by simulating J/ψ → e + e − decays, passing them through GEANT3 simulation [38] of the STAR detector, embedding the simulated digital signals into real data, and finally reconstructing the embedded events through the same procedure as for the real data. The central values and statistical errors of the J/ψ polarization parameters are obtained by maximizing the likelihood and corrected for possible biases that are estimated from a toy Monte Carlo (ToyMC). In this ToyMC, the same numbers of J/ψ signal and background candidates as in real data are randomly generated with fixed values of polarization parameters after applying detector acceptance and reconstruction efficiencies. The extracted J/ψ polarization parameters from the pseudo-data following the same procedure as described above are compared to the input values in terms of both central values and statistical errors, and the differences are applied as corrections to real data, which are generally very small compared to statistical errors.

C. Signal extraction
Invariant mass spectra of electron-positron pairs are shown in Fig. 1 for five different p J/ψ T bins. The combinatorial background contribution is estimated from samesign charge pairs, shown as filled areas in the figure. The raw numbers of J/ψ candidates are estimated by subtracting same-sign distributions from opposite-sign ones and integrating resulting counts within the invariant mass window of 3 − 3.15 GeV/c 2 . The contribution from the residual background is found to be between 1.5-2.5% and thus neglected here. The two-dimensional N J/ψ (cos θ, φ) distributions in the HX and CS frames are shown in Fig. 2. The J/ψ reconstruction efficiency multiplied by the detector acceptance, A × ε(cos θ, φ), are shown in Fig. 3, corresponding to the invariant mass window of 3 − 3.15 GeV/c 2 . The detector acceptance, track reconstruction, BEMC electron identification and HT trigger efficiencies are estimated from simulation. Polarization of input J/ψ's do not play any role due to two-dimensional determination of the efficiencies. The electron identification efficiencies due to application of TPC and TOF requirements are estimated from data [22] using a pure electron sample from gamma conversions. The electron dE/dx and 1/β distributions are fit with a Gaussian distribution to calculate the cut efficiencies. The TOF matching efficiency is evaluated using TPC tracks that are matched to BEMC hits in order to suppress the pileup contribution. The bias due to the geometrical correlation between BEMC and TOF acceptance is corrected using an electron sample from data.
To check the results obtained from fit, the uncorrected J/ψ distributions are compared to the expected ones as shown in Figs. 4 and 5. The former are obtained by projecting two-dimensional N J/ψ (cos θ, φ) distributions onto either the cos θ or φ direction, while the latter are generated using the extracted J/ψ polarization parameters from data and taking into account the detector acceptance and efficiency. The expected J/ψ distributions agree well with the measured ones, confirming that the maximum likelihood method can be used to reliably extract the J/ψ polarization parameters. Also shown in these figures as references are the expected cos θ and φ distributions corresponding to the extreme cases where the polarization parameters (λ θ , λ φ , λ θφ ) = (±1, 0, 0) and (λ θ , λ φ , λ θφ ) = (0, ±1, 0) are used, respectively.
The systematic uncertainties are estimated for the following sources. 1) Acceptance: in extracting efficiencies from simulation, different parameterizations of the inclusive J/ψ p T and rapidity spectra [22] are tried and the difference is used as the uncertainty.
2) PID: the uncertainty in the electron identification GeV/c). The black markers (blue filled histograms) are the spectra from opposite-sign (same-sign) charge pairs, while the red markers represent those obtained by subtracting the same-sign spectra from the opposite-sign ones. efficiencies is assessed by varying the mean and width of the TPC nσ e and TOF 1/β distributions according to their uncertainties in the efficiency calculation, and by simultaneously varying the cut in both data and simulation on the ratio between the track momentum and energy deposition in the BEMC from 0.3 < pc/E < 1.5 to 0.2 < pc/E < 1.4 or 0.4 < pc/E < 1.6. Additional uncertainties are considered in evaluating the TOF matching efficiency, including a correction factor to account for the correlation between BEMC and TOF acceptances in obtaining the TOF matching efficiency from data.
3) Tracking: the uncertainty in track reconstruction efficiency is obtained by simultaneously varying the cuts in data and simulation on the minimum number of TPC hit points from 20 to 18, 19, 22, or 25, on maximum DCA from 1 cm to 0.8 or 1.2 cm, and by varying the momentum resolution in simulation within its uncertainty. 4) Triggering: the uncertainty in the HT trigger efficiency is obtained by simultaneously changing the HT trigger threshold cut with ±5% variation.
For each of the systematic sources, the same analysis procedure is followed and the resulting maximum differences to the default results are taken as the uncertainties. The total systematic uncertainties are a quadrature sum of individual sources, as shown in Table I. IV. J/ψ → µ + µ −

A. Dataset, event and track selections
The dataset was taken for p+p collisions at √ s = 200 GeV in 2015, and corresponds to an integrated luminosity of 122 pb −1 . Events are selected online with a dimuon trigger, which requires at least two signals in the MTD whose timing difference to the start time provided by the VPD falls within the pre-defined trigger timing window.
Events used in offline analysis are required to have a vertex position of |V TPC z | < 100 cm along the beam direction to maximize statistics. To reject the events originating from the beam hitting the beam pipe, primary vertices are further required to be within 2 cm radially with respect to the center of the beam pipe.
In the analysis of the dimuon decay channel, charged tracks reconstructed in the TPC should have at least 15 TPC space points used for reconstruction. The ratio of the actually used to the maximum possible number of TPC space points is required to be larger than 0.52 to reject split tracks. The distance of closest approach (DCA) to the primary vertex needs to be less than 3 cm to suppress contribution from secondary decays and pileup tracks. The selected TPC tracks are afterwards refit with the primary vertex included in order to improve the momentum resolution. Tracks are then propagated from the TPC to the MTD radius. Only tracks with p T > 1.3 GeV/c are selected to achieve high efficiency for reaching the MTD after losing energy along the trajectory. Once a track is matched to the closest MTD hit, cuts on variables, ∆y, ∆z and ∆t tof , are applied to further suppress background hadrons. Here, ∆y and ∆z are the residuals between the projected track position at the MTD radius and the matched MTD hit along azimuthal and beam directions, respectively. We require ∆y and ∆z to be within 3 (3.5) σ of their resolutions for p T < (>) 3 GeV/c. ∆t tof is the difference between the measured time-of-flight with the MTD and the calculated time-offlight from track extrapolation with a muon particle hypothesis, and should satisfy |∆t tof | < 1 ns. Additional PID capabilities arise from the energy loss measurement in the TPC. It is quantified as nσ π , whose definition is similar to that of electrons as described in Sect. III A but using a pion hypothesis. A cut of −2 < nσ π < 3 is applied.

B. Analysis procedure
To extract the J/ψ polarization in the dimuon decay channel, a different strategy is adopted compared to the one used for the dielectron channel as described in Sect. III B. Equation 1 is integrated over φ and cos θ, yielding two 1-D distributions: and  The λ θφ term vanishes in both integrations. The polarization parameters, λ θ and λ φ , are extracted from a simultaneous fit to corrected J/ψ yield distributions as a function of cos θ and φ of daughter µ + with Eqs. 5 and 6. This strategy is motivated by the worse signal-to-background ratio for the dimuon decay channel compared to the dielectron decay channel, and fitting the 1-D distributions of Eqs. 5 and 6 is therefore more stable. However, the λ θφ parameter cannot be extracted from this method.  The number of J/ψ extracted in each cos θ or φ bin needs to be corrected for the detector acceptance and efficiency, denoted as A × ε(cos θ, φ). It is evaluated via simulation as described in Sect. III B but for the muon channel. Since A × ε(cos θ, φ) depends on both cos θ and φ, the projected 1-dimentional(1-D) A × ε as a function of cos θ or φ is affected by the assumed polarization of input J/ψ in the simulation. On the other hand, the λ θφ value does not affect the averaged A×ε as A×ε(cos θ, φ) is symmetric with respect to cos θ = 0 and |φ − π/2| = 0. Given that the J/ψ polarization is not known a priori and the correction for the detector acceptance and efficiency depends on it, an iterative procedure is adopted. In the first iteration, the 1-D A × ε as a function of cos θ or φ is evaluated using non-polarized J/ψ in the simulation, and the polarization parameters are extracted from data after correcting for A × ε. In the second iteration, the extracted polarization parameters from the previous iteration are used in the simulation to assess A × ε, which in turn is used to correct data and obtain new polarization parameters. The iteration continues until the differences of the obtained polarization parameters between two con-secutive iterations are less than 0.01. This threshold is determined based on the statistical precision of the data.
To validate the iterative procedure, a ToyMC is developed which is different from that used in the electron channel analysis. The single muon efficiency as a function of p T , η and φ, extracted from the GEANT simulations, is applied to mimic realistic detector acceptance and detection efficiency. J/ψ's with realistic λ θ and λ φ values in four different p T bins, as presented in Sec. V, are used as input to the ToyMC while the λ θφ value is assumed to be 0. Both pseudo-data and J/ψ pseudo-efficiency are generated in the ToyMC. Depending on the statistical precision of the pseudo-data and how the pseudo-efficiency is obtained, the following tests are done: • Test 1 -Large statistics with correct efficiency ("Large stat., corr. eff."): the pseudo-data sample has significantly larger statistics than real data, and the pseudo-efficiency is generated using the same polarization parameters as for pseudo-data. This represents a best-case scenario, and the polarization parameters, λ θ , λ φ and λ inv , are extracted using Eqs. 5 and 6. Differences to the input polarization values are shown in Fig. 6 as open circles for both HX and CS frames. In most cases, the input values are recovered with small discrepancies arising from the limited acceptance of the MTD.
• Test 2 -Limited statistics with correct efficiency ("Limited stat., corr. eff."): the pseudo-data sample has comparable statistical precision to real data, while the pseudo-efficiency is generated using the same polarization parameters as for pseudo-data.
To avoid random fluctuation of one pseudo-data sample, 500 independent samples of similar statistics are generated. The mean values of the polarization parameters extracted from the 500 pseudodata samples are compared to the input values, and the differences are shown in Fig. 6 as filled squares. Compared to "Test 1", the extracted polarization values deviate further from the input ones due to the influence of the limited statistics in the pseudodata sample on top of the limited MTD acceptance.
• Test 3 -Limited statistics using the iterative procedure ("Limited stat., last iteration"): in the last test, the 500 pseudo-data samples are generated with comparable statistical precision to real data, but the iterative procedure as described above is used to obtain the efficiency. Polarization values equal to 0 are used in the first iteration, and the procedure stops after the same convergence criteria of 0.01, as for the real data, between two consecutive iterations is fulfilled. The resulting differences to the input values are shown as open squares in Fig. 6, and agree with "Test 2" quite well. This indicates that very small biases are introduced in the iterative procedure.
It has also been found that the correct polarization values are always obtained as long as the convergence occurs and no matter what input polarization values are used in the first iteration. The ToyMC validation confirms that the J/ψ polarization parameters can be extracted reliably using the iterative procedure. The residual biases shown in Fig. 6 are corrected for, as described in Sect. IV C.

C. Signal extraction
The selected muon candidates of opposite-sign charges are paired, and the resulting invariant mass distribution is shown in Fig. 7 for the entire sample used in the dimuon channel analysis. The raw J/ψ yield is extracted by fitting the invariant mass distribution with a Gaus-sian function describing the J/ψ signal, and a polynomial function describing the background. Data points in the ψ(2S) mass region (3.6 < M µµ < 3.8 GeV/c 2 ) are excluded from the fit. The mean of the Gaussian distribution is fixed to the J/ψ mass in the PDG [39]. In total, the J/ψ yields are extracted in ten cos θ bins and fifteen φ bins for each p J/ψ T interval. The order of the background polynomial function ranges from 2 to 5, depending on the p J/ψ T and the µ + cos θ and φ bin. To facilitate the fits below 2 GeV/c, the widths of the Gaussian function in individual cos θ and φ bins are fixed to be the same as that extracted from fitting the inclusive invariant mass distribution integrated over cos θ and φ bins in the same p J/ψ T interval. For p J/ψ T above 2 GeV/c, the width of the Gaussian function is left as a free parameter. Variations in the following aspects of the fit procedure are applied: the bin width of the invariant mass distribution, fixing the width of the Gaussian function also for p J/ψ T above 2 GeV/c, the order of the polynomial function and the fit range. The average J/ψ yields from these variations are used for extracting the polarization parameters. J/ψ yields with significance less than 3 are not considered. Upper panels of Fig. 8 show an example of the average raw J/ψ yield, depicted as open circles, as a function of cos θ and φ for 0 < p J/ψ T < 1 GeV/c in the HX frame.
Following the iterative procedure, the efficiency multiplied by the detector acceptance from the last iteration is shown in the upper panel of Fig. 8 as dashed lines. It is scaled to the same integral as the data distribution. The lower panels of Fig. 8 show the fully corrected J/ψ yield as a function of cos θ and φ, along with the simultaneous fit to both distributions using Eqs. 5 and 6 (solid lines). The polarization parameters, λ θ and λ φ , are obtained from the simultaneous fit and listed in the figure. Similar plots in the CS frame are shown in Fig.  9 for 0 < p J/ψ T < 1 GeV/c. The average λ θ and λ φ values from 24 combinations of different track quality and muon identification cuts are taken as the central values. Figures 8 and 9 show results from one such combination as an example. As shown in Fig. 6, small biases are present in the extracted polarization parameters, due to a combination of limited MTD acceptance, limited statistical precision of data and the usage of the iterative procedure. Using the ToyMC described above, the extracted J/ψ polarization parameters are compared with input values in terms of central values and statistical errors, and the differences are applied as corrections to real data, which are also very small compared to statistical errors.
The systematic uncertainties arise from signal extraction, detector acceptance and efficiency. As mentioned above, different aspects of the fitting procedure are varied, and 24 combinations of track quality and muon identification cuts are tried. For the latter, cuts are changed consistently in data and MC simulation, and the entire analysis chain is repeated for each case. The RMS's of θ cos  the distributions for polarization parameters stemming from these two groups of variations are taken as the systematic uncertainties, respectively. The total systematic uncertainties are a quadrature sum of individual sources, as shown in Table II.

V. RESULTS
The polarization parameters, λ θ , λ φ and λ θφ for inclusive J/ψ are measured in both the HX and CS frames in p+p collisions at √ s = 200 GeV/c, as shown in Fig. 10. The dimuon and dielectron results are shown as filled and open symbols respectively, and are consistent with each other in the overlapping p T range even though they cover different rapidity regions. All three polarization parameters are consistent with 0 within statistical and systematic uncertainties, except for λ θ in the CS frame above 8 GeV/c whose central value is at -0.69 ± 0.22 ± 0.07. No strong p T dependence is seen in all cases. The numerical values of the measured J/ψ polarization parameters are listed in Appendix (Tables IV−VII). Model calculations for prompt J/ψ from ICEM [7], NRQCD with two sets of LDMEs denoted as "NLO NRQCD1" [40] and "NLO NRQCD2" [13], are shown in Fig. 10 for comparison. The prompt J/ψ's include directly produced J/ψ and those from decays of excited charmonium states. Also shown in Fig. 10 are CGC+NRQCD calculations for direct J/ψ. The feeddown contributions from excited charmonium states and b-hadron decays, present in data, are expected to have a small impact on the theoretical calculation of the J/ψ polarization parameters [19]. For λ θ in the HX frame, the ICEM calculation predicts a sizable transverse polarization at low p T , while the J/ψ polarization from CGC+NRQCD changes from slightly transverse at low p T to slightly longitudinal at higher p T . The difference between the LDMEs used in the two NLO NRQCD calculation is that additional η c production data measured by the LHCb Collaboration [14] is used to determine LDMEs for "NLO NRQCD1" besides those used for the case of "NLO NRQCD2". They show opposite behaviors for λ θ and λ φ in both reference frames. To quantify the agreement between data and model calculations, the χ 2 test has been performed simultaneously using the data points in HX and CS frames for both channels. The χ 2 /NDF and corresponding p-values are listed in Table  III. While no model can be ruled out definitively based solely on the data presented, the CGC+NRQCD gives the best overall description. The λ inv values extracted according to Eq. 2 for inclusive J/ψ are shown in Fig. 11 as a function of p T for both the HX and CS frames. The dimuon and dielectron results are shown as filled and open circles, respectively. The vertical bars represent the statistical errors while the boxes around data points depict the systematic uncertainties. The λ inv values measured in the two frames are consistent with each other within experimental uncertainties, confirming the reliability of the results. The λ inv values are consistent with the CGC+NRQCD calculations within uncertainties.

VI. SUMMARY
For the first time, the inclusive J/ψ polarization parameters, λ θ , λ φ and λ θφ , are measured as a function of p T in p+p collisions at √ s = 200 GeV in both the Helicity and Collins-Soper reference frames. Results utilizing the dimuon and dielectron decay channels are presented and agree with each other within uncertainties although slightly different kinematic ranges are covered. The inclusive J/ψ's do not exhibit significant transverse or longitudinal polarization with little dependence on p T .
Among several model calculations compared to data, the CGC+NRQCD agrees the best overall. These results provide additional tests and valuable guidance for theoretical efforts towards a complete understanding of the J/ψ production mechanism in vacuum. frames. Open and filled symbols are for measurements through the dielectron and dimuon decay channels covering different rapidity ranges. The vertical bars represent the statistical errors while the boxes around data points depict the systematic uncertainties. The data points are placed at the mean pT value determined from the inclusive p J/ψ T spectrum measured in p+p collisions at √ s = 200 GeV [22]. Model calculations [7,13,19,40] are also shown for comparison. The ICEM and two NLO NRQCD calculations are for prompt J/ψ, while the CGC+NRQCD is for direct J/ψ.  11. λinv of J/ψ vs. pT in both HX (circles) and CS (diamonds) reference frames. The open and filled symbols are for measurements through the dielectron and dimuon decay channels, respectively. The vertical bars represent the statistical errors while the boxes around data points depict the systematic uncertainties. CGC+NRQCD [19] calculations are also shown for comparison. IV. The inclusive J/ψ polarization parameters in the HX frame in different pT bins measured through the dielectron channel within |y| < 1. The first uncertainty is statistical and the second is systematic.