Estimates for the single-spin asymmetries in $p^{\uparrow}p \to J/\psi X$ process at PHENIX RHIC and SPD NICA

We study the transverse single-spin asymmetry (TSSA) in $p^{\uparrow}p \to J/\psi X$ reaction, incorporating both transverse-momentum and spin effects. To predict production cross section of prompt $J/\psi$ we use two different approaches, the non-relativistic QCD (NRQCD) factorization approach and the Improved Color Evaporation Model (ICEM), and show how the predicted results for TSSAs depend on choice of hadronization model. For initial-state factorization we consider two models: the standard Generalized Parton Model (GPM) and the Colour Gauge-Invariant version of it (CGI-GPM). We demonstrate that PHENIX collaboration data on TSSA in the process $p^{\uparrow}p \to J/\psi X$ constrain the gluon Sivers function of the proton and rule-out one of existing parameterizations. Estimates for the TSSAs in $p^{\uparrow}p \to J/\psi X$ process for the conditions of the future SPD NICA experiment are presented for the first time.


I. INTRODUCTION
The transverse-momentum dependent (TMD) parton distribution functions (PDFs) incorporate information about three-dimensional structure of proton and it's spin properties [1,2]. Among the leading-twist TMD PDFs, the Sivers function [3,4] is one of the most interesting and it is widely investigated in p ↑ p → hX inclusive processes [5,6]. Sivers function describes the number density of unpolarized gluons g (or quarks q) with intrinsic transverse-momentum q T inside a transversely polarized proton p ↑ , with three-momentum P and spin polarization vector S, where x is the proton light-cone momentum fraction carried by the gluon, F g (x, q T ) is the unpolarised TMD parton density, ∆ N F ↑ g (x, q T ) is the Sivers function, q T = |q T | and symbol (ˆ) denotes a unit vector,â = a/|a|.
In the present paper we are interested in accessing the gluon Sivers function using the TSSA of inclusive J/ψ-production. However, the task of studying gluon distributions using charmonia is rather challenging theoretically. Production of charmonia proceeds in two stages: first, a cc-pair is produced at short distances, predominantly via gluon-gluon fusion but also with a non-negligible contribution of qq−and qg−initiated subprocesses. The second stage is hadronization of cc-pair into a physical charmonium state, which proceeds essentially nonperturbatively, at large distances (low scales) and is accompanied by a complicated rearrangement of color via exchanges of soft gluons between the cc-pair and other colored partons produced in the collision. At present, two approaches to describe cc hadronization are most popular: the non-relativistic QCD (NRQCD) factorization [7] and (Improved) Color-Evaporation Model (CEM) [8][9][10].
In a context of TMD-factorization, first rigorous results for heavy quarkonium physics have been obtained only very recently [11,12], showing that the TMD-factorization formula for quarkonium production will differ from the case of Drell-Yan pair or Higgs-boson production processes. For quarkonia, the TMD-factorization formula have to include additional shape-functions with the corresponding evolution. However, in the present paper we adopt a simpler phenomenological approach of Generalized Parton Model (GPM) which will be described in more details in Sec. II A. In the standard GPM approach the Sivers function is assumed to be process-dependent, because effects of initial (ISI) and final-state interactions (FSI) are factorized into this function. An alternative approach is Color-Gauge-Invariant (CGI) GPM formalism of Refs. [16-19, 21, 22], which we summarize in Sec. II B and II C. In this framework the process-dependent ISIs and FSIs are lifted from Sivers-like TMD PDF to the coefficient function, using the one-gluon exchange approxiamtion. Thus in CGI-GPM results for Sivers function extracted from different processes can be directly compared.
The behavior of TSSA in the process p ↑ p → J/ψX has been studied recently in Refs. [19,21,23]. In Refs. [19,21], J/ψ production was treated in color-singlet approximation [19] and full NRQCD approach [21], including color-octet states, respectively. In case of colorsinglet mechanism, the CGI-GPM corrections to GPM cross sections have been included. In Ref. [23], the ICEM was used and authors investigated the effect of the evolution of TMD PDFs involved, on the asymmetry A N .
The aim of the present study, is to to provide reasonable estimates for TSSA-effects which could be observed in the kinematic conditions of planned Spin Physics Detector experiment at NICA collider [13,14]. We will compare estimates for TSSA, obtained in GPM and CGI-GPM frameworks and will explore both theoretical approaches for charmonium hadronization -NRQCD and ICEM. In case of NRQCD, we will show that at small transverse momenta of produced charmonium k T C < m C , color-singlet mechanism describes data for prompt J/ψ production and inclusion of color-octet terms is not needed at LO in α s . Therefore for our NRQCD predictions for TSSA within GPM and CGI-GPM we include only color-singlet states of final cc-pair. In connection to this findings, we comment on the results of of Refs. [19,21] in Sec. III.
The present paper is organized as follows. In the Sec. II, we present basic formulas of GPM and CGI-GPM for 2 → 1 and 2 → 2 partonic subprocesses, as well as details on NRQCD approach and ICEM. In the Sec. III, our numerical results for transverse-momentum spectra of prompt J/ψ and TSSA A N in p ↑ p → J/ψX at PHENIX RHIC and SPD NICA are presented and discussed. In this section we also compare our results with the similar results obtained earlier for PHENIX RHIC at √ s = 200 GeV both in NRQCD approach [19,21] and ICEM [23].

A. GPM and CGI-GPM
The factorization scheme, which is suitable to describe inclusive observables in hadronic collisions, depends on hierarchy between typical hard scale of the process and involved transverse momenta. Collinear Parton Model (CPM) is used for studies of high-p T production, when hard scale (µ ∼ p T Λ QCD ). Such a way, in CPM we neglect transverse momenta of partons initiating the hard process and hadronic cross section can be factorized as convolution of hard coefficients and collinear parton distribution functions (PDFs) f a (x, µ F ) at the factorization scale µ F with µ F µ ∼ p T and a = q,q, g. In the kinematical domain of CPM, transverse momenta of final-state particles (k T ) are generated in the hard scattering of partons and influence of small intrinsic transverse momenta of partons in hadrons (q T ), which originates from nonperturbative effects, can be neglected. The typical estimate for average squared intrinsic transverse-momentum of a parton is q 2 T 1 GeV 2 or even smaller.
If one is interested in particle production with small transverse momenta k T q 2 T µ in a hard processes with the scale µ, effects of intrinsic parton motion in hadrons have to be taken into account. In the Transverse Momentum Dependent (TMD) factorization, the hadronic cross section is expressed as a convolution of TMD PDFs F a (x, q T , µ, ζ) and hard partonic cross section. The TMD factorization theorem has been proven in the limit of k T µ or q T µ, and TMD PDFs evolve with respect to two scales: µ and ζ, where the latter one is the so-called rapidity scale related with rapidity divergences [24]. The typical value of hard scale µ in charmonium (C) production is given by charmonium mass, GeV, so TMD factorization should be used in the region k T m C . The region k T m C requires matching with finite-order perturbative corrections and possible nonperturbative power-suppressed corrections to the TMD term.
The Generalized Parton Model (GPM) is a simplified version of TMD factorization, which is generally applied for phenomenological estimates of various observables in proceses for which TMD factorization have not been rigorously proven yet. Typically TMD PDFs in GPM are parametrized by a simple factorized prescription: where f a (x, µ F ) is corresponding collinear PDF. The dependence on transverse-momentum of a parton is described by Gaussian distribution G a (q T ) = exp[−q 2 T / q 2 T a ]/(π q 2 T a ) with normalization condition d 2 q T G a (q T ) = 1. And effects of TMD-evolution w.r.t. rapidity scale [24,25] ζ are neglected, which means that GPM estimates are applicable only in a narrow range of hard and rapidity scales ζ ∼ µ. The latter condition is however always fulfilled in our case, since the scale of the process of charmonium production at low k T C is given by m C .
Within the GPM, the differential cross section for charmonium production in protonproton collisions for the 2 → 1-type hard subprocess g(q 1 ) + g(q 2 ) → C(k) can be written as where C = J/ψ, ψ(2S), or χ c (1P ), and with s = 2P 1 P 2 -the squared center of mass energy of the pp-collision. Integrating-out delta functions one obtains: For consistency of GPM and to avoid problems with gauge invariance of hard scattering amplitudes, four-momenta of initial-state partons have to be put on mass-shell (q 2 1 = q 2 2 = 0). Hence q 1,2 read: where x 1,2 are the proton light-cone momentum fractions carried by the partons: As follows from Eqs. (6) and (7): The latter result allows one to integrate-out delta function δ(ŝ − m 2 C ) taking the integral over dx 2 in Eq. (5). To perform the calculation of transverse-momentum spectrum in fixed rapidity interval one inserts a Heaviside theta-function implementing the rapidity cut under the sign of integral (5), while the rapidity of C can be calculated as The cross section for the 2 → 2 subprocess (g(q 1 )+g(q 2 ) → C(k)+g(q 3 ), C = J/ψ, ψ(2S)) is given by formula (3) with Replacing q 2 3 →ŝ +t +û − m 2 C witht = (q 1 − k) 2 ,û = (q 2 − k) 2 one can remove integral over d 4 q 3 by delta function and obtain Then one can take integral over dx 2 analytically using delta function δ(ŝ +t +û − m 2 C ). In such a way, no kinematic approximations is made and exact 2 → 1 and 2 → 2 kinematics is implemented in presence of transverse-momentum of initial-state partons.
Master formulas presented above have been used directly in calculations of prompt J/ψ production in the NRQCD approach, see Sec. (II B). In the case of ICEM approach the treatment of 2 → 1 subprocess is slightly different, see Sec. II C.
In this paper we study TSSAs, usually denoted by A N , measured in p ↑ p → CX (C = J/ψ, χ c , ψ(2S)) inclusive reactions and defined as: where ↑, ↓ are opposite proton spin orientations perpendicular to the scattering plane in pp center-of-mass frame. The numerator and denominator of A N reads whereF ↑,↓ g (x, q T , µ F ) is the distribution of unpolarized gluon (or quark) in polarized proton. Following the Trento conventions [26], the gluon Sivers function(GSF) can be introduced as (14) and GSF has to satisfy the positivity bound where M p -mass of the proton.
We adopt factorized Gaussian parametrizations for both the unpolarized TMD distribution F g (x, q T , µ F ) and the Sivers function which satisfies the bound (15) for any values of α and β. After introducing the parameter we write for gluon Sivers function (GSF): GSF set  In our numerical calculations we will use two different of GSF obtained earlier in Refs. [6] which we call SIDIS1 for brevity, and [27] which we refer to as GSF parametrization by Corresponding values of parameters are collected in the Table I.
To introduce the CGI-GPM let us first recall the explanation of Sivers effect, which had been described for the first time in Ref. [20]. In this paper it was shown, that Sivers asymmetry in semi-inclusive DIS (SDIS) process at leading twist is a quantum effect generated by exchanges of soft gluons between initial (ISI) or final-state (FSI) partons produced in a hard process, and spectator system originating as a remnant of an incoming hadron. In standard TMD factorization, this soft gluons are taken into account within the gauge-invariant definition of Sivers-like TMD PDF, which contains Wilson lines. The sign of Sivers TMD PDF depends on the direction of Wilson lines, which can be past-or future-pointing, representing the "space-time trajectory" of an initial-state or struck quark produced respectively in Drell-Yan or SDIS hard-scattering process. Thus the Sivers function in standard TMD and perhaps in GPM approaches is process-dependent and it is not clear how to extend factorization for Sivers effect to the processes with colored final-states, like J/ψ production.
The aim of CGI-GPM [16-19, 21, 22] formalism is to extract above-mentioned processdependence from the TMD PDF to the hard-scattering coefficient. The effects of ISI and FSI are included in CGI-GPM via one-gluon exchange approximation [16,19]. For the case of gluon Sievers effect, this approximation leads to appearance of independent GSFs of f - 1T ) corresponding to two independent ways of combining three gluons into a color-singlet (Fig. 1). The coupling of additional "eikonal" gluon from the GSF to the hard process leads only to modification of the color structure of the latter one.
There is no four-momentum transfer from the additional gluon to the hard process, because Sivers effect comes from imaginary part of the loop integral over momentum of exchanged gluon [20], while the latter one is saturated by the contribution of the soft region. Moreover, only coupling of eikonal gluon to the initial-state or observed colored final-state particles contributes to the asymmetry, while effects of coupling to un-observed final-state partons cancel-out between amplitude and complex-conjugate amplitude. While arguments above are specific for the one-gluon exchange approximation, an additional argument in favor of CGI-GPM is, that it's hard-scattering coefficients reproduce coefficients in the twist-3 collinear approach [16]. In the Fig. 1 we collect the corresponding Feynman rules and prescriptions for calculation of the hard scattering coefficient for the numerator (12) of the TSSA within CGI GPM. For detailed derivation see e.g., Ref. [19]. The color projectors for f and d-type GSFs are defined as follows: where f abc (d abc ) is totally antisymmetric (symmetric) structure constant of the SU (N c ) color gauge group, t a ij -the generators of SU (N c ) group in the fundamental representation, and with N c = 3.
In the following sections we mostly use the CGI-GPM hard-scattering coefficients which already had been calculated by other authors, in which case we cite the corresponding reference. However in Sec. II C we will need new (to our knowledge) CGI-GPM hardscattering coefficient for gg → cc-process with both final-state c-quarks being observed.

B. NRQCD
In the framework of the NRQCD-factorization approach, the cross section of charmonium production via a partonic subprocess a+b → C +X is given by a double expansion in powers of α s and squared relative velocity of heavy quarks in a bound state v 2 as: where n denotes the set of color, spin, orbital and total angular momentum quantum numbers of the cc pair and the four-momentum of the latter is assumed to be equal to the one of the physical quarkonium state C. The cross section dσ(a + b → cc[n] + X) can be calculated in perturbative QCD as an expansion in α s . The nonperturbative transition of the cc pair LDMEs can be determined from measured decay widths of charmonia using the known nextto-leading-order (NLO) QCD result or from calculations in potential models [28]. Coloroctet LDMEs are considered as free parameters in charmonium production cross sections.
Typically LDMEs up to NNLO (O(v 4 )) in v 2 -scaling are included in NRQCD-factorization where J = 0, 1, 2. For the LDMEs of P -vawe states we adopt the known Heavy-Quark Spin Symmetry relations, which are valid up to O(v 2 ): From general point of view, LDMEs should be universal and process-independent parameters. However in practice their numerical values strongly depend on approach which is used to describe cc−pair production and data included into the fit. For example one can compare color-octet LDMEs obtained in LO CPM [29,30], in NLO CPM [31] and in the k T −factorization approach [32,33]. Neevertheless, the hierarchy expected from velocity-scaling rules: ] is respected by all the fits. We adopt the following values of color-singlet LDMEs [34]: Squared LO in α S amplitudes for 2 → 1 subprocesses in CPM are well-known [29]: As it was discussed above, in GPM these subprocesses contribute to transversemomentum spectrum because TMD PDFs are involved.
There is only one relevant LO in α S 2 → 2 partonic subprocess, which describes direct 1 ] as it is in the Color-Singlet Model, because 2 → 2 subprocesses producing color-octet states of cc-pair are of NLO in α s in the GPM approach for k T C -spectrum. The squared amplitude for this partonic subprocess reads [35]: Turning now to the case of CGI-GPM factorization we introduce the following notations: for the coefficient-function, obtained within the CGI-GPM factorization prescription. Then, following Ref. [19], one writes-down the contribution of the 2 → 1 or 2 → 2 subprocess of production of 3 P 1 -states of cc-pair to the numerator of TSSA as: where F g(f ) 1T and F g(d) 1T are above-mentioned f -type (C-even) and d-type (C-odd) GSFs, and ⊗ denotes convolution in the light-cone momentum fraction and transverse-momentum of gluon from the polarized proton. Here C U is the color factor of the unpolarized cross section, which corresponds to the usual QCD result, while C  [16-19, 21, 22]: for the case of 3 S 1 final-state (2 → 2 process), while states (2 → 1 processes). One notices, that in both cases, d-type GSF does not contribute, so that only f -type GSF is relevant for color-singlet model.  A comment on the treatment of color-octet contributions in the CGI-GPM approach of Refs. [21,22] is in order here. As it is explained in Ref. [22] around Eq. (A34) the coefficient functions contributing to the numerator of the asymmetry formula for 2 → 1subprocesses producing color-octet states of cc-pair in this approach are equal to zero due to cancellation between ISI and FSI contributions, unlike coefficient functions of ordinary GPM in our Eqs. (29), (30) and (32). Due to this fact, authors of Refs. [21,22] had to include contributions of 2 → 2-processes to obtain non-zero effect on the asymmetry from color-octet channels. However, since LO in α s contribution in CGI-GPM is zero, one expects that NLO 2 → 2 contribution should not contain any infra-red or collinear divergences. Contrary to this expectation, e.g. coefficient function for the contribution of g + g → cc 1 S clearly contains non-integrable singularities att → 0 orŝ → m 2 C which due to non-zero transverse-momentum of initial-state partons will lead to divergent cross section in GPM.
This fact has been mentioned in Ref. [21], where it is admitted that a regulator µ IR ∼ 0.8 GeV for this divergences have to be introduced. In our opinion, appearance of non-regulated divergences in CGI-GPM deserves further study and we are going to address this problem in the future. In the numerical calculations of the present paper we include only color-singlet NRQCD channels, which are free from above-mentioned problem. We should also point-out, that usage of color-octed LDMEs of NRQCD, obtained in the NLO fits together with LO coefficient functions is not consistent, because NLO corrections to short-distance coefficients in NRQCD are very significant and even change sign of P -wave contributions. Nevertheless, such NLO LDMEs had been used in Refs. [21,22]. The good agreement of un-polarized cross section calculated in full NRQCD with NLO LDMEs in Ref. [21] with experimental data is probably due to neglecting of the feed-down contributions by these authors, while this contributions are actually non-negligible (see our Tab. II). Also, the treatment of 2 → 1 kinematics in Refs. [21,22] is different from ours, which leads to significant numerical effects for NICA energies.

C. ICEM
Main physical assumption of the ICEM is that all cc−pairs with invariant masses below the DD-threshold hadronize to charmonia with some probability, which is independent from angular momentum and spin quantum numbers of the cc-pair. In the ICEM [9, 10] the invariant mass of the intermediate charm quark-antiquark pair is constrained to be larger than the mass of produced charmonium state, m C , instead of using the same lower limit of integration -2m c as it was done in the traditional CEM [8]. As a result, the ICEM describes the charmonium yields as well as the ratio of ψ(2S) over J/ψ better than the old CEM.
The partonic cross section, differential in cc-invariant mass is related to the well-known total cross section of production of cc-pairs as a function of partonic squared center-of-mass energyŝ: with w = 4m 2 c /ŝ, as follows: so that the GPM-factorization formula for production of cc-pairs with invariant mass M and total three-momentum k via gluon-gluon fusion can be written as: where the invariantŝ = k 2 = (q 1 + q 2 ) 2 can be represented as in Eq. (8). Finally, for differential cross section of chamonium C production in proton-proton collision in ICEM one has: where F C is the process-independent hadronization probability to the charmonium state C.
Then one integrates-out q 2T and M 2 using delta functions to find In this equation the integral over dx 2 also can be removed by delta function δ(q 3 1 +q 3 2 −k 3 ) thus obtaining the master formula for numerical calculations. The quark-antiquark annihilation channel for cc−pair production has been incorporated into the calculation in a similar way.
In case of CGI-GPM factorization, the numerator of the TSSA (12), reads: whereσ (f /d) CGI (ŝ, gg → cc) is the f /d-type coefficient function of CGI-GPM integrated over phase-space of final-state cc-pair with fixed invariant massŝ as in Eq. (38). Since in the ICEM both c-andc-quarks in the final state are observed, the corresponding hard-scattering coefficient is different from the coefficient for D-meson TSSA, given e.g. in Ref. [37]. To obtain new coefficient functions we take into account interactions of eikonal gluon with the initial-state gluon coming from un-polarized proton as well as with final-statec and c-quarks (middle and right diagrams of the Fig. 2). The f /d-type hard-scattering coefficients thus obtained have the form: Integrating this coefficient-functions over the phase-space of the final-state with fixed cc invariant-massŝ one obtains: It is interesting, that integrated hard-scattering coefficient for d-type GSF is equal to zero similarly to the case of NRQCD, so that in both of our models, heavy-quarkonium TSSA is sensitive only to f -type GSF.
To obtain prompt-J/ψ production spectra we take into account direct as well as feed-down contributions from decays of χ cJ and ψ(2S)-states. At the stage of numerical calculation in ICEM, we put m c = 1.2 GeV and charmonium masses are taken from PDG tables: m J/ψ = 3.096 GeV, m ψ(2S) = 3.686 GeV, m χ c0 = 3.415 GeV, m χ c1 = 3.510 GeV, and m χ c2 = 3.556 GeV.

A. PHENIX RHIC
To begin with, we compare theoretical predictions obtained in NRQCD-factorization approach with recent experimental data for transverse-momentum spectra of prompt J/ψmesons, measured by PHENIX RHIC experiment [38]. In our NRQCD calculations we take the charm quark mass as one half of mass of the physical charmonium state, m c = m C /2, while in the ICEM calculations it is kept fixed at m c = 1.2 GeV. Also, in case of feed-down production, the kinematic effect of the mass splittings between charmonium states turns out to be significant and we take into account momentum-shift between high-mass charmonium state and final J/ψ meson as it was done e.g. in Ref. [39]: k T J/ψ = (m J/ψ /m C )k T C .

Phenomenological analysis of intrinsic transverse-momentum of partons in proton in LO
and NLO of CPM [40] demonstrates that for gluon one has q 2 T 1 GeV 2 and the same estimation was obtained for J/ψ production in GPM [19].
Throughout our analysis the renormalization and factorization scales has been identified and chosen to be µ F = µ R = ξ k 2 T + m 2 C where ξ is varied between ξ = 1/2 and ξ = 2 about its default value ξ = 1 to estimate the theoretical uncertainty due to the freedom in the choice of scales. The resulting errors are indicated as shaded bands in our figures, however they mostly cancel-out in asymmetries A N .
The direct production of J/ψ-mesons at the O(v 0 ) includes only contributions from CSM (34). The color-octet states 3 P  does not contribute if initial-state partons are on-mass shell. Further color-octet contributions (29), (28), and (33) are suppressed as O(v 4 ) so it is natural to expect them to be negligible at small k T C < m C . In fact, similarly to the results of Ref. [19], we have found that taking into account only color-singlet production mechanism, the good description of prompt J/ψ transverse-momentum spectra at small k T J/ψ < m J/ψ , i.e. in the region of applicability of TMD-factorization, can be achieved in GPM see the left panel of Fig. 3.
Also, our NRQCD calculation leads to total cross section ratios of direct and feed-down contributions in good agreement with experimental data of Ref. [38], see Tab. II.   However, we should emphasize, that our calculations are different from calculations of the Ref. [19] in a respect, that we consistently take into account feed-down contributions from ψ(2S) and χ cJ -states, while in this reference they where added very crudely, by multiplication Very similar predictions for transverse-momentum spectrum can be obtained in the framework of ICEM, see the right panel of Fig. 3. The values of hadronization probabilities used are F J/ψ = 0.02, F χ c1 = F χ c2 = 0.06 and F ψ(2S) = 0.08. They had been obtained via the fit of total cross section of J/ψ-production at PHENIX and above-mentioned experimentallymeasured fractions of J/ψ-feeddown contribution form ψ(2S) and χ cJ -decays (Tab. II). This values of hadronisation probabilities are numerically close to the values obtained in Ref. [10] by the fit of LHC data in the k T -factorization approach for cc-pair production.
Our estimations for TSSAs at PHENIX kinematic conditions, obtained in the GPM accompanied by NRQCD-factorization approach or ICEM, are shown by thin histograms in the Fig. 4 and Figs. 5-6 as functions of x F and transverse-momentum respectively, together with the recent experimental data from Ref. [41]. We conclude that within standard GPM initial-state factorization, the parametrisation for Sivers function by D'Alesio et al. is marginally consistent with experimental data for both hadronization models, while SDIS1parametrisation predicts too large effects at positive x F ψ and is essentially ruled-out for the case of ordinary GPM initial-state factorization.
The TSSA results for the CGI-GPM initial-state factorization are presented in the same

B. SPD NICA
In this section we present our predictions for J/ψ transverse-momentum spectrum and TSSA in the kinematic conditions of planned SPD NICA experiment in proton-proton collisions with √ s = 24 GeV. The SPD is expected to be an almost 4π-geometry detector [13][14][15], thus a relatively wide coverage in rapidity |y| < 3 can be achieved.
As for k T J/ψ -spectrum, the GPM calculations, both in NRQCD-factorization and ICEM, lead to results consistent with NRQCD predictions of Parton Reggeization Approach [42] at small transverse-momentum, while the latter predictions are in agreement with NLO turn out to be consistent with PHENIX data. Thus we conclude that we can safely perform predictions for TSSA at NICA energies.
Estimates for TSSA at SPD NICA experiment, computed within NRQCD and ICEM approaches under GPM initial-state factorization assumption are shown in Figs. 8 and 9 respectively for the SDIS1 and D'Alesio et al. parametrisations for GSF. We find that for standard GPM initial-state factorization, the SDIS1 predicts gigantic values for asymmetries at NICA energies -up to 60% (Fig. 8). However such big effects can hardly be expected to appear, since this parametrisation contradicts PHENIX data, when GPM is used. GPM predictions with D'Alesio et al. parametrisation (Fig. 9) look more realistic and they are quite robust against the choice of J/ψ-formation model. Measurable asymmetries up to 5% for the x F ψ -spectrum and up to 2% for k T J/ψ -spectrum are predicted.
Our results obtained using CGI-GPM initial-state factorization are shown in Figs

Conclusions
In the present paper we have performed a phenomenological analysis of gluon Sivers function contribution to the transverse TSSA of prompt J/ψ-production within NRQCDfactorization (essentially Color-Singlet Model in our case) and ICEM for the description of J/ψ-formation, employing both state-of-art initial-state factorization models: GPM and CGI-GPM. The goal of our analysis was to make predictions for TSSA in the kinematic conditions of planned SPD NICA experiment. We have found, that within standard GPM initial-state factorization, the SDIS1 parametrisation for gluon Sivers function contradicts PHENIX data, while parametrisation of D'Alesio et al. leads to reasonable predictions for magnitude, J/ψ transverse-momentum and x F ψ dependence of the asymmetry with |A N | < ∼ 2 − 3% (Fig. 9). Within CGI-GPM initial-state factorization, contradiction of SDIS1-parametrization with PHENIX data is eliminated, and it predicts |A N | < ∼ 5 − 10% at SPD NICA kinematics (Fig. 10). Hence, observation of sizable transverse TSSA in in- clusive J/ψ-production does not contradict existing experimental data and their theoretical interpretation within a wide range of J/ψ-formation and initial-state factorization models.
In any case, measurements at SPD NICA will significantly constrain our knowledge about gluon Sivers function in a proton.