Single-spin asymmetry of J/ψ production in p+p, p+Al, and p+Au collisions with transversely polarized proton beams

C. Aidala, Y. Akiba, 51, ∗ M. Alfred, V. Andrieux, N. Apadula, H. Asano, 50 B. Azmoun, V. Babintsev, A. Bagoly, N.S. Bandara, K.N. Barish, S. Bathe, 51 A. Bazilevsky, M. Beaumier, R. Belmont, A. Berdnikov, Y. Berdnikov, D.S. Blau, 42 M. Boer, J.S. Bok, M.L. Brooks, J. Bryslawskyj, 8 V. Bumazhnov, S. Campbell, V. Canoa Roman, R. Cervantes, C.Y. Chi, M. Chiu, I.J. Choi, J.B. Choi, † Z. Citron, M. Connors, 51 N. Cronin, M. Csanád, T. Csörgő, 62 T.W. Danley, M.S. Daugherity, G. David, 56 K. DeBlasio, K. Dehmelt, A. Denisov, A. Deshpande, 56 E.J. Desmond, A. Dion, D. Dixit, J.H. Do, A. Drees, K.A. Drees, J.M. Durham, A. Durum, A. Enokizono, 52 H. En’yo, S. Esumi, B. Fadem, W. Fan, N. Feege, D.E. Fields, M. Finger, M. Finger, Jr., S.L. Fokin, J.E. Frantz, A. Franz, A.D. Frawley, Y. Fukuda, C. Gal, P. Gallus, P. Garg, 56 H. Ge, F. Giordano, Y. Goto, 51 N. Grau, S.V. Greene, M. Grosse Perdekamp, T. Gunji, H. Guragain, T. Hachiya, 51 J.S. Haggerty, K.I. Hahn, H. Hamagaki, H.F. Hamilton, S.Y. Han, J. Hanks, S. Hasegawa, T.O.S. Haseler, X. He, T.K. Hemmick, J.C. Hill, K. Hill, A. Hodges, R.S. Hollis, K. Homma, B. Hong, T. Hoshino, N. Hotvedt, J. Huang, S. Huang, K. Imai, M. Inaba, A. Iordanova, D. Isenhower, D. Ivanishchev, B.V. Jacak, M. Jezghani, Z. Ji, X. Jiang, B.M. Johnson, 20 D. Jouan, D.S. Jumper, J.H. Kang, D. Kapukchyan, S. Karthas, D. Kawall, A.V. Kazantsev, V. Khachatryan, A. Khanzadeev, C. Kim, 30 E.-J. Kim, M. Kim, D. Kincses, E. Kistenev, J. Klatsky, P. Kline, T. Koblesky, D. Kotov, 53 S. Kudo, K. Kurita, Y. Kwon, J.G. Lajoie, A. Lebedev, S. Lee, S.H. Lee, 56 M.J. Leitch, Y.H. Leung, N.A. Lewis, X. Li, S.H. Lim, 63 M.X. Liu, V-R Loggins, S. Lökös, 17 K. Lovasz, D. Lynch, T. Majoros, Y.I. Makdisi, M. Makek, V.I. Manko, E. Mannel, M. McCumber, P.L. McGaughey, D. McGlinchey, 34 C. McKinney, M. Mendoza, A.C. Mignerey, D.E. Mihalik, A. Milov, D.K. Mishra, J.T. Mitchell, G. Mitsuka, S. Miyasaka, 58 S. Mizuno, 59 P. Montuenga, T. Moon, D.P. Morrison, S.I. Morrow, T. Murakami, 50 J. Murata, 52 K. Nagai, K. Nagashima, T. Nagashima, J.L. Nagle, M.I. Nagy, I. Nakagawa, 51 K. Nakano, 58 C. Nattrass, T. Niida, R. Nouicer, 51 T. Novák, 62 N. Novitzky, A.S. Nyanin, E. O’Brien, C.A. Ogilvie, J.D. Orjuela Koop, J.D. Osborn, A. Oskarsson, G.J. Ottino, K. Ozawa, 59 V. Pantuev, V. Papavassiliou, J.S. Park, S. Park, 54, 56 S.F. Pate, M. Patel, W. Peng, D.V. Perepelitsa, 12 G.D.N. Perera, D.Yu. Peressounko, C.E. PerezLara, J. Perry, R. Petti, M. Phipps, 24 C. Pinkenburg, R.P. Pisani, M.L. Purschke, P.V. Radzevich, K.F. Read, 57 D. Reynolds, V. Riabov, 49 Y. Riabov, 53 D. Richford, T. Rinn, S.D. Rolnick, M. Rosati, Z. Rowan, J. Runchey, A.S. Safonov, T. Sakaguchi, H. Sako, V. Samsonov, 49 M. Sarsour, S. Sato, B. Schaefer, B.K. Schmoll, K. Sedgwick, R. Seidl, 51 A. Sen, 57 R. Seto, A. Sexton, D. Sharma, I. Shein, T.-A. Shibata, 58 K. Shigaki, M. Shimomura, 41 T. Shioya, P. Shukla, A. Sickles, C.L. Silva, D. Silvermyr, B.K. Singh, C.P. Singh, V. Singh, M.J. Skoby, M. Slunečka, M. Snowball, R.A. Soltz, W.E. Sondheim, S.P. Sorensen, I.V. Sourikova, P.W. Stankus, S.P. Stoll, T. Sugitate, A. Sukhanov, T. Sumita, J. Sun, Z Sun, Z. Sun, J. Sziklai, K. Tanida, 51, 54 M.J. Tannenbaum, S. Tarafdar, 61 A. Taranenko, A. Taranenko, 55 G. Tarnai, R. Tieulent, 36 A. Timilsina, T. Todoroki, M. Tomášek, C.L. Towell, R.S. Towell, I. Tserruya, Y. Ueda, B. Ujvari, H.W. van Hecke, J. Velkovska, M. Virius, V. Vrba, 26 N. Vukman, X.R. Wang, 51 Y.S. Watanabe, C.P. Wong, C.L. Woody, C. Xu, Q. Xu, L. Xue, S. Yalcin, Y.L. Yamaguchi, 56 H. Yamamoto, A. Yanovich, J.H. Yoo, I. Yoon, H. Yu, 48 I.E. Yushmanov, W.A. Zajc, A. Zelenski, S. Zharko, and L. Zou

We report the transverse single-spin asymmetries of J/ψ production at forward and backward rapidity, 1.2 < |y| < 2.2, as a function of J/ψ transverse momentum (pT ) and Feynman-x (xF ). The data analyzed were recorded by the PHENIX experiment at the Relativistic Heavy Ion Collider in 2015 from p+p, p+Al, and p+Au collisions with transversely polarized proton beams at √ s N N = 200 GeV. At this collision energy, single-spin asymmetries for heavy-flavor particle production of p+p collisions provide access to the spin-dependent gluon distribution and higher-twist correlation functions inside the nucleon, such as the gluon Qiu-Sterman and trigluon correlation functions. Proton+nucleus collisions offer an excellent opportunity to study nuclear effects on the correlation functions. The data indicate a positive asymmetry at the two-standard-deviation level in the p+p data for 2 GeV/c < pT < 10 GeV/c at backward rapidity, and negative asymmetries at the twostandard-deviation level in the p+Au data for pT < 2 GeV/c at both forward and backward rapidity, while in p+Al collisions the asymmetries are consistent with zero within the range of experimental uncertainties.

I. INTRODUCTION
In polarized p+p collisions, the transverse single spin asymmetry (TSSA), A N , is defined as the amplitude of the azimuthal angular modulation of the outgoing particle's scattering cross section with respect to the transverse spin direction of the polarized proton. Early theoretical predictions which were purely based on perturbative calculations showed that the TSSA should be inversely proportional to the hard scale of the scattering [1], and if applied to the reactions at energies accessible at the Relativistic Heavy Ion Collider (RHIC), the asymmetry would be very small, of order 10 −4 . However, the experimental asymmetries of light-flavored hadrons turned out to be much larger, A N = O(10 −1 ) [2,3].
To explain what has been observed in experiments, several theoretical frameworks [4][5][6][7][8] were developed in the 1990s. In the collinear factorization framework, contributions from multi-parton correlations to the transverse-spin-dependent cross section were introduced through three types of spin-momentum correlations: (1) twist-3 correlation functions of a polarized hadron φ (3) a/A (x 1 , x 2 , s ⊥ ) convolved with leading-twist correlation functions of an unpolarized hadron φ b/B (x ) and with leading-twist parton fragmentation functions D c→C (z), (2) transversity parton distribution functions δq a/A (x, s ⊥ ) convolved with twist-3 correlation functions of an unpolarized hadron φ (3) b/B (x 1 , x 2 ) and leadingtwist parton fragmentation functions D c→C (z), and (3) transversity parton distribution functions δq a/A (x, s ⊥ ) convolved with leading-twist correlation functions of an unpolarized hadron φ b/B (x ) and with twist-3 fragmentation functions D (3) c→C (z 1 , z 2 ) [9]: (1) In this notation, a/A means the distribution of parton a in hadron A, b/B means the distribution of parton b in hadron B and c → C means the fragmentation of parton c into hadron C. Additionally, x is the Bjorken * PHENIX Spokesperson: akiba@rcf.rhic.bnl.gov † Deceased parton momentum fraction of the incoming hadron; z is the fraction of the outgoing partonic momentum carried by the detected hadron; s ⊥ is the transverse spin of the incoming hadron;σ,σ , andσ are the partonic hardscattering cross sections of a process where higher twist is associated with the incoming polarized or unpolarized hadron, or outgoing partons, respectively. Non-zero values for φ a/A (x 1 , x 2 , s ⊥ ) corresponds to initial-state effects [8] and D (3) c→C (z 1 , z 2 ) to final-state effects [10]. Initial-state effects are described by the twist-3 three-parton correlation functions φ which measure the quantum interference between two scattering amplitudes of the incoming hadron [11], while final-state effects are related to twist-3 fragmentation functions D (3) c→C (z 1 , z 2 ) which describe the process in which the outgoing parton fragments into a final-state hadron [11]. At RHIC energies, heavy quark production, such as J/ψ production, is dominated by gluon-gluon interactions. Because the gluon transversity distribution does not exist, the second and third terms of Eq. 1 are zero. This means that heavy flavor A N is free from finalstate effects and is sensitive to initial-state effects, such as the gluon Qiu-Sterman and trigluon correlations which correspond to the factor φ (3) a/A (x 1 , x 2 , s ⊥ ) in Eq. 1 [12]. In the case of high energy hadronic collisions, a nonvanishing SSA is generated by a parton-level spin flip and a phase difference between the scattering amplitude and the corresponding complex conjugate. In the context of the quantum chromodynamic (QCD) collinearfactorization framework, the parton-level spin flip is generated by the interference between the active single parton and a two-parton composite state of the scattering amplitude. On the other hand, the phase difference is achieved by the interference between the real and imaginary part of the partonic scattering amplitude [8,9,13]. For the Qiu-Sterman correlation, the quantum interference is between a quark state of momentum fraction x and a quark-gluon composite state with the same momentum fraction where either the gluon or quark carries the total momentum of the quark-gluon composite state [14]. In the trigluon correlation, the two parton composite state is composed of two gluons instead of a quark and gluon as described above for the Qiu-Sterman correlation [11,15]. The collinear factorization framework has been widely used to describe the TSSA's measured at RHIC [16][17][18][19][20][21][22][23].
An alternative treatment is known as the Transverse-Momentum-Dependent (TMD) formalism. In this for-malism, the cross section is factorized into hardscattering cross sections and TMD parton distribution and fragmentation functions (PDFs and FFs) [24]. For the TMD approach to be valid in the context of p+p collisions, Q 2 must be large, in order to use perturbative QCD, while the transverse momentum must satisfy p T Q and not be much larger than the intrinsic, parton transverse momentum k T , so that effects of the latter remain visible [25]. One of the TMD PDFs, called the Sivers function [4], is widely used in describing the TSSA's that were observed in different processes [26][27][28][29][30]. The Sivers function, denoted by f ⊥ 1T (x, k 2 ⊥ ), describes the distortion in the distribution of unpolarized partons with momentum fraction x and transverse momentum k ⊥ in a transversely polarized hadron. This distorted distribution of unpolarized partons causes an azimuthal anisotropy in the distribution of parton transverse momenta in the polarized hadron which gives rise to the nonzero TSSA. As it has been described above, at low p T , the nonperturbative TMD Sivers function will be responsible for its SSA, while twist-3 dominates the contributions to the SSA when p T ∼ Q. At intermediate p T , one can see the transition between these two frameworks and a relation between Sivers Function and Qiu-Sterman Function has been shown in Ref. [31].
With increasing experimental information on the quark Sivers function during the last ten years, our understanding of this quantity has matured [32][33][34][35][36], while the gluon Sivers function is still relatively unknown. The transversely polarized p+p collisions studied at RHIC present a very good opportunity to study the gluon Sivers function as gluon-gluon interactions are dominant in p+p collisions at RHIC energies. PHENIX has measured the TSSA for J/ψ production at central and forward rapidities [37] and, at small p T values, where the J/ψ mass becomes the large scale Q in TMD factorization, the result has been compared to a gluon Sivers function derived in the context of the color-evaporation model in [38] and the generalized-parton model in [39].
In proton-nucleus (p+A) collisions, the increase of the atomic number results in increasing gluon occupancy and therefore gluon saturation effects may become important in the small x region. In the Regge-Gribov limit, the properties of saturated gluons in the infinite-momentum frame can be described by the Color Glass Condensate (CGC) which has been applied to a variety of processes such as e+p, e+A, A+A, and p+A collisions [40]. The quark and gluon distribution functions for large nuclei were computed first in Ref. [41] in the weak coupling limit. Using the CGC framework, one can describe the rescatterings of the outgoing parton within the nucleus. In the coherent QCD multiple scattering framework [42], it has been shown that at low p T the rates of single and double hadron production are highly suppressed and the amount of suppression grows with rapidity and centrality; meanwhile, at high p T , such nuclear modification effects become less pronounced as we enter the perturbative region in a dilute nuclear medium. For the com-putation of TSSA's, a hybrid approach has been widely used in photon, γ-jet and dijet production [43][44][45][46]. The hybrid approach treats the gluon distribution inside the heavy nucleus in the CGC framework and utilizes the twist-3 formalism for the proton side.The resummation of the power corrections in p+A reactions shifts the nuclear PDF to higher x [42]. The new PHENIX p+A collision data offer the opportunity to quantify how this shift in x affects the twist-3 description of the TSSA's discussed above.

II. EXPERIMENT SETUP
The J/ψ mesons are measured by detecting muon pairs µ + µ − in the two PHENIX muon spectrometers. The two muon spectrometers with full azimuthal coverage cover a range in pseudorapidity η ∈ [−2.2, −1.2] for the south arm and η ∈ [1.2, 2.4] for the north arm (see Fig. 1). The momenta of the muons are measured by the muon trackers (MuTr). The MuTr is composed of three stations of cathode-strip tracking chambers inside a radial magnetic field [47]. The MuTr chambers are followed by the muon identifiers (MuID) which contain 5 sensitive layers (named gap0-gap4) per arm, with each layer made of one vertically and one horizontally oriented Iarocci tube plane, all interleaved with 10-or 20-cm thick steel absorber plates to suppress hadron backgrounds. With the front absorber provided by the PHENIX central magnets and the MuID together with additional stainless steel absorber mounted in 2011, the total thickness of the steel absorbers is about 190 cm. The probability for a 4 GeV/c pion to punch through the whole absorber is less than 3% [48,49]. Muons with momenta of at least 3 GeV/c can penetrate all the absorbers and reach the last layer of the MuID with high efficiency.
Two beam-beam counters (BBC) are located at opposite ends at 144 cm from the interaction point along the beam line with full azimuthal coverage and |η| ∈ [3.1, 3.9]. Each BBC has 64 elements consisting of quartzČerenkov radiator and mesh-dynode PMT. The BBCs also serve as an interaction vertex finder with resolution along the beam direction of about 2 cm and the z-vertex approximately follows a Gaussian distribution centered at 0 with a width of about 10 cm in p+p collisions and in addition play the role of a luminosity detector [50]. For the minimum-bias (MB) trigger, it is required to have one or more hits in each BBC. The MB-trigger efficiency for p+p collisions is about 55% [51], while this trigger is 84% (72%) efficient for p+Au (p+Al) collisions; the percentage is defined with multiplicity in the south BBC (Agoing direction) [52].
Events containing dimuon candidates were selected using the combination of BBC-MB trigger and other muon triggers. The "2-deep muon trigger" requires that both muon tracks have at least one hit in the last two MuID gaps and no less than two hits in other gaps. The "sagitta-3 muon trigger" selects tracks that were recorded by MuTr with high momentum, requiring that the maximum track sagitta, determined by the 3 MuTr stations, be less than 3 MuTr cathode strips at the middle plane.
The dimuon candidates are then screened by the dimuon cuts shown in Table I. Here, p z and p T are respectively the longitudinal and transverse momentum of the dimuon pair with respect to the beam direction, the lower-side p T cut is imposed due to the MuTr resolution, DG0 is the distance between two matched tracks coming one each from the MuTr and MuID, projected to the first MuID gap, DDG0 is the opening angle between these projected tracks, ntrhits (nidhits) is the number of hits generated by a track in the MuTr (MuID), lastgap is the last-gap number in the MuID that a track penetrates, and DCA r is the distance of closest approach between a muon track and the beam line. The event vertex determined by the BBC follows a Gaussian distribution centered at 0 with a width of about 40 cm. A z-position vertex cut within ±30 cm ensures that the collision occurs inside the spectrometer acceptance and keeps approximately 50% of all collisions. The vertex χ 2 is determined by fitting the two candidate tracks with the event vertex given by the BBC.
In the 2015 RHIC run, we took 11 weeks of p+p collision data with an integrated luminosity of about 40 pb −1 and with both the blue beam (clockwise) and yellow beam (counter-clockwise) transversely (vertically) polarized at the PHENIX interaction point. In p+A collisions, only the blue beam (the proton beam) was polarized; we have 2-and 5-week data sets for p+Al and p+Au with integrated luminosities of about 6.0 pb −1 and 6.6 pb −1 , respectively. Table II shows the beam polarizations. Each proton beam was filled into 111 bunches with different spin states in 8 base patterns [53].

III. DATA ANALYSIS
The maximum likelihood method was used to extract the transverse single spin asymmetry A N . The likelihood L for one dimuon pair with azimuthal angle φ with re-spect to the incoming polarized proton beam direction is 1 + P · A N sin(φ pol − φ), where P is the beam polarization, and φ pol is the direction of beam polarization which is +(−)π/2 when the spin is up (down). Then the log-likelihood for all the dimuon pairs is given by: We select the value of A N that maximizes L. The statistical uncertainty of A N is obtained by calculating the inverse of the second derivative of L with respect to A N : Fit of inclusive A(φ) with a cosine modulation. The asymmetry and χ 2 per degree of freedom from the fit are shown.
2.8 GeV/c 2 and 3.4 GeV/c 2 . Then we corrected that A N with the background A N (A bgr N ) using Eq. 4 and calculated the corresponding uncertainty in A J/ψ N using Eq. 5: Here, f is the background fraction in the the J/ψ invariant mass region. The statistical uncertainty from f is not taken into account because its contribution to the total statistical uncertainty is small (< 5%). The background A N under the J/ψ peak was estimated from a sideband in the dimuon invariant mass range of 1.5-2.4 GeV/c 2 ; a similar approach was used in previous PHENIX measurements [54,55]. The background fraction f was obtained using the Gaussian Process Regression (GPR) method [56][57][58][59] for each p T and x F bin. The GPR method is a nonparametric regression approach considered to be less biased, and has been used in several previous PHENIX measurements [55,60]. Figure 2 shows the fitting result with the GPR method using p+Au data. The dimuon invariant mass windows of 1.5 GeV/c 2 to 2.2 GeV/c 2 and 4.3 GeV/c 2 to 6.0 GeV/c 2 are used for estimating the background shape (see Fig. 2). Then the signal (J/ψ and ψ(2S)) yields are obtained by subtracting the GPR background from the inclusive dimuon spectrum. The J/ψ and ψ(2S) peaks are fitted respectively by a Crystal Ball function [61] and a Gaussian function. The vertical lines in Fig. 2 represent ±2σ mass windows around the J/ψ peak, where the ψ(2S) contribution will be negligible; J/ψ production is dominant in both the p-going and Au-going directions. The background fraction is 17% for the Au-going direction, 13% for the Al-going direction and 10% for the p-going direction.
The A incl N and A bgr N asymmetries for each p T and x F bin are calculated with the maximum likelihood method described above but with different dimuon invariant mass windows. For A incl N we used unlike-sign muon pairs in the invariant mass range ±2σ around the J/ψ peak, while for A bgr N we used the fixed invariant mass range from 1.5 GeV/c 2 to 2.4 GeV/c 2 .
In this analysis, there are two sources of systematic uncertainty, detailed as sources 1 and 2 in the following two paragraphs and listed quantitatively in Tables III  and IV. For the total systematic uncertainty displayed in Figs. 4 and 5, we have combined these two sources quadratically.
The first systematic uncertainty source ("source 1") concerns the method of determining the asymmetry itself. We check this by determining A N with a different method, the azimuthal fitting method. Similar to Eq. 2, the production cross section of J/ψ as a function of the azimuthal angle φ is given by: where φ pol = +(−)π/2 when the spin is up (down). The asymmetry can be written as function of azimuthal angle φ as: Therefore, A N can be extracted by fitting the A(φ) with a cosine modulation. As an example, Figure 3 shows the determination of A incl N for dimuons with 0.42 < p T < 2 GeV/c in the Au-going direction. The differences of J/ψ A N determined from the maximum likelihood method and azimuthal fitting method are treated as a systematic uncertainty. The value of the source 1 uncertainty ranges from 1% to 35% of the statistical uncertainties.
A second source of systematic uncertainty ("source 2") is from the method of determining the background fraction f . We studied a potential bias of the GPR method by parameterizing the background with a 3rd order polynomial instead. The f changed by about 2% and the corresponding difference in the resulting background corrected J/ψ A N s has been assigned as a systematic uncertainty which is of the order of 10% of the statistical uncertainty. Table III show [37] within one-standarddeviation. The 2015 p+p data favor a positive asymmetry (at the 2σ level) in the high-p T bin in backward rapidity. With limited statistics, the A N in all p T and x F bins for p+Al collisions are consistent with zero. In p+Au collisions, the asymmetry in the high-p T bin is consistent with zero, although there is a trend to a nonzero A N (at the 2σ level) in the low-p T bin in both the forward and backward directions. Figure 5 and Table IV show the A J/ψ N as a function of x F . In the p+p data, it is consistent with the previous PHENIX results [37] and a ∼ 2σ positive A N is observed in the backward higher x F bin. The result for the other x F bins are consistent with zero. For the p+Au data, a ∼ 2σ negative A N is observed in the forward high-x F bin and the backward low-x F bin. A scale uncertainty from the polarization (3%) is not included in both Figure 4 and Figure 5.

V. CONCLUSION
We have reported the measurements of the transverse single-spin asymmetry in J/ψ production at forward and backward rapidity with 1.2 < |y| < 2.