Probing $J/\psi$ Production Mechanisms in Proton–Proton Collisions at SPD/NICA Energies

Abstract

We investigate inclusive $J/\psi$ production in proton–proton collisions at tens of GeV $\sqrt{s}$ energy, relevant for forthcoming measurements with the Spin Physics Detector (SPD) at NICA. Simulations are performed using the PEGASUS event generator with transverse-momentum-dependent (TMD) gluon densities, comparing the recent KMR-based KL${}^{\prime}2025$ and CCFM-based LLM${}^{\prime}2024$ parametrizations. Differential cross sections in rapidity and transverse momentum exhibit smooth, stable behavior under renormalization-scale variation, while factorization-scale dependence exposes limitations of the LLM${}^{\prime}2024$ set at low scales in contrast to KL${}^{\prime}2025$. Normalized $p_{T}$ spectra reveal distinct hardening patterns linked to the underlying gluon $k_{T}$ broadening in each model. The relative contributions of color-singlet and color-octet channels are also quantified, demonstrating the dominance of color-octet mechanisms in the SPD energy regime. These results provide the first detailed assessment of quarkonium production sensitivity to gluon TMDs near threshold, offering timely theoretical guidance for upcoming $J/\psi$ measurements at SPD/NICA.

1. Motivation— Quarkonium production in hadronic collisions serves as a powerful probe of Quantum Chromodynamics (QCD) across its perturbative and nonperturbative domains ALICE:2021dtt; Arbuzov:2020cqg. In particular, the production of the charmonium state $J/\psi$ is directly sensitive to gluon dynamics: at leading order it proceeds via gluon–gluon fusion, while its hadronization into a bound state involves the interplay between short-distance scattering and long-distance QCD effects Arbuzov:2020cqg; physics5030044; Petrelli:1997ge.

Measurements of inclusive $J/\psi$ production over a broad range of center-of-mass energies have enabled stringent tests of theoretical frameworks such as the Color-Singlet Model (CSM) Baier:1983va, Nonrelativistic QCD (NRQCD) factorization Bodwin:1994jh; Butenschoen:2012px, and the Color-Evaporation Model (CEM) Frawley:2008kk. High-precision LHC measurements ATLAS:2011aqv; LHCb:2021pyk; ALICE:2021dtt have established a detailed picture of quarkonium production at multi-TeV energies, where global analyses combining color-singlet (CS) and color-octet (CO) mechanisms achieve good agreement with data.

At moderate collision energies, however, the situation is less constrained. Here, the relevant gluon momentum fractions are larger, the available phase space for high-$p_{T}$ recoils is reduced, and the sensitivity to the transverse motion of partons becomes clearer Aidala:2012mv. In this regime, the transverse-momentum-dependent (TMD) formalism provides a natural extension of QCD factorization, allowing direct access to the intrinsic $k_{T}$ structure and polarization of gluons inside the proton physics5030044; Arbuzov:2020cqg. Quarkonium production, particularly inclusive $J/\psi$, thus emerges as a promising observable to probe gluon TMDs and their evolution.

The forthcoming Spin Physics Detector (SPD) experiment at the NICA collider will explore proton–proton collisions at $\sqrt{s}\leq 27~\mathrm{GeV}$ with polarized beams Arbuzov:2020cqg; physics5030044. This energy range bridges the gap between fixed-target and high-energy collider experiments, providing a uniquely clean environment to investigate gluon dynamics in the transition region where perturbative and nonperturbative effects overlap SPDproto:2021hnm.

In this work, we investigate inclusive $J/\psi$ production in proton–proton collisions for $\sqrt{s}\leq 27~\mathrm{GeV}$ using the PEGASUS event generator, which implements TMD gluon densities in a $k_{T}$-factorized framework Lipatov:2019oxs. Two recent TMD parameterizations are employed: the KMR-based KL${}^{\prime}2025$ Kotikov:2025wft and the CCFM-based LLM${}^{\prime}2024$ Lipatov:2024xni sets. We present differential cross sections in transverse momentum and rapidity, including theoretical uncertainties from renormalization-scale variations. The normalized $p_{T}$ spectra reveal distinct behaviors for the two TMD sets, reflecting their different evolution patterns at SPD scale. Additionally, the relative contributions of CS and CO channels are analyzed across beam energies, elucidating how production mechanisms evolve from the SPD regime to the high-energy domain. These results establish the predictive foundations for quarkonium measurements at SPD/NICA and quantify the sensitivity of $J/\psi$ observables to the underlying gluon densities.

2. Theoretical Framework— Within the NRQCD factorization approach, the inclusive $J/\psi$ production cross section is written as a sum over intermediate $c\bar{c}$ states $n$  Bodwin:1994jh; Petrelli:1997ge,

$$d\sigma(pp\to J/\psi+X)=\sum_{n}d\hat{\sigma}(pp\to c\bar{c}[n]+X)\,\langle\mathcal{O}^{J/\psi}[n]\rangle,$$ (1)

where $d\hat{\sigma}$ represents the perturbatively calculable short-distance coefficient (SDC) for the production of a heavy-quark pair in state $[n]$, and $\langle\mathcal{O}^{J/\psi}[n]\rangle$ denotes the long-distance matrix element (LDME) governing its nonperturbative transition into a physical $J/\psi$ meson.

At leading order in $\alpha_{s}$, two classes of partonic subprocesses contribute:

1. CS mechanism:

   $$g^{*}+g^{*}\to c\bar{c}\big[{}^{3}S_{1}^{(1)}\big]+g,$$ (2)

   where the $c\bar{c}$ pair is produced directly in a CS configuration with the same quantum numbers as the $J/\psi$ Petrelli:1997ge. The final-state gluon ensures color conservation and carries part of the recoil momentum.

2. CO mechanisms:

   $$g^{*}+g^{*}\to c\bar{c}[n],\quad n\in\big\{{}^{1}S_{0}^{(8)},\,{}^{3}S_{1}^{(8)},\,{}^{3}P_{J}^{(8)}\,(J=0,1,2)\big\},$$ (3)

   where the $c\bar{c}$ pair is produced in a color-octet state and subsequently evolves into the physical $J/\psi$ through soft-gluon emissions, modeled via electric dipole (E1) transitions Petrelli:1997ge. These $2\to 1$ subprocesses contribute at the same order in $\alpha_{s}$ as the CS $2\to 2$ process.

Both CS and CO channels are implemented in the PEGASUS Lipatov:2019oxs event generator within the TMD factorization framework, where the initial gluons are off-shell ($g^{*}$) and carry intrinsic transverse momentum $k_{T}$. The generated event samples include all relevant intermediate $c\bar{c}[n]$ states ($[{}^{3}S_{1}^{(1)}]$, $[{}^{1}S_{0}^{(8)}]$, $[{}^{3}S_{1}^{(8)}]$, $[{}^{3}P_{J}^{(8)}]$), allowing a consistent treatment of both color and spin dynamics of $J/\psi$ production at SPD/NICA energies.

In this study, the LDMEs are set to unity to isolate the kinematic dependence of the SDCs. The partonic cross sections are convoluted with TMD gluon densities evolved via the CCFM and KMR schemes. The renormalization and factorization scales are chosen as $\mu^{2}=m_{T}^{2}=M_{J/\psi}^{2}+p_{T}^{2}$, with the $J/\psi$ mass consistent with PDG values ParticleDataGroup:2018ovx; Lipatov:2019oxs.

Monte Carlo event generation and phase-space integration are performed using PEGASUS, which provides .lhe samples containing weighted events according to the $k_{T}$-factorization formalism Lipatov:2019oxs; Alwall:2006yp. Twelve statistically independent runs, each with $4\times 10^{5}$ events, are produced to ensure stable predictions of differential distributions in $p_{T}$, $y$, and $\sqrt{s}$. For consistent comparison across different TMD gluon sets, the results are normalized to the total hadronic cross section where needed.

This framework provides a unified and quantitative basis to explore the sensitivity of inclusive $J/\psi$ production to the underlying gluon dynamics in the SPD/NICA energy range.

Figure 1: Rapidity distributions ${d\sigma}/{dy}$ for inclusive $J/\psi$ production at $\sqrt{s}=9$, $18$, and $27~\mathrm{GeV}$ obtained using the KL${}^{\prime}2025$ Kotikov:2025wft and LLM${}^{\prime}2024$ Lipatov:2024xni gluon densities. The shaded bands represent the uncertainty from renormalization-scale variation. A clear broadening and enhancement of the distributions with increasing $\sqrt{s}$ reflect the expanding rapidity phase space and growing gluon–gluon luminosity at higher energies.

Figure 2: Differential cross sections ${d\sigma}/{dp_{T}}$ for inclusive $J/\psi$ production in pp collisions at $\sqrt{s}=9$, $18$, and $27~\mathrm{GeV}$ using the KL${}^{\prime}2025$ Kotikov:2025wft and LLM${}^{\prime}2024$ Lipatov:2024xni gluon densities. The shaded bands indicate the uncertainty due to renormalization-scale variation. The comparison illustrates the evolution of the $p_{T}$ spectrum with increasing collision energy.

3. Results and Discussion— Figure 1 shows the rapidity distributions $d\sigma/dy$ for inclusive $J/\psi$ production at $\sqrt{s}=9$, $18$, and $27~\mathrm{GeV}$. The spectra exhibit a symmetric and centrally peaked structure across all studied energies, with the peak position remaining close to mid-rapidity ($|y_{\text{peak}}|<0.07$), consistent with the expected symmetry of proton–proton collisions in this regime.

A systematic broadening of the rapidity distribution with increasing energy is observed for both KL${}^{\prime}2025$ and LLM${}^{\prime}2024$ sets, as reflected by the rising RMS and FWHM values. The width roughly doubles from $\sqrt{s}=9$ to $27~\text{GeV}$, indicating the progressive opening of the rapidity phase space and the expanding kinematic reach in gluon momentum fractions. Concurrently, the overall magnitude of the distribution increases, reflecting the enhanced gluon–gluon luminosity that drives quarkonium production from the near-threshold to the intermediate-energy regime accessible at SPD.

The comparison between the two gluon density parameterizations reveals that LLM${}^{\prime}2024$ systematically gives higher cross sections at all rapidity values, while maintaining a slightly broader shape than KL${}^{\prime}2025$. This difference originates from the broader intrinsic transverse-momentum width and the softer small-$x$ gluon behavior encoded in the LLM${}^{\prime}2024$ evolution, which enhance central production rates. At the lowest energy, $\sqrt{s}=9~\mathrm{GeV}$, both models show a mild suppression around $y\simeq 0$, indicative of near-threshold kinematic constraints on the accessible gluon flux.

The shaded bands represent the theoretical uncertainty estimated by varying the renormalization scale around its central value, $\mu_{R}=m_{T}$, by a factor of two, i.e., $\mu_{R}=m_{T}/2$ and $\mu_{R}=2m_{T}$, while keeping the factorization scale fixed at $\mu_{F}=m_{T}$. The quantitative spread and overall behavior of these bands remain stable across the considered energy range. It is worth noting that variation of the factorization scale does not serve as an equally reliable estimator of theoretical uncertainty in this regime, as the LLM${}^{\prime}2024$ set exhibits a rapid suppression of the cross section under downward scale shifts, indicating the onset of a kinematic threshold in its low-scale applicability.

The two gluon-density sets thus yield nearly identical rapidity shapes but notably different normalizations, indicating that the overall production rate is primarily controlled by the gluon-luminosity normalization, whereas the shape is determined by the kinematic mapping of the partonic subprocess. This systematic comparison highlights that, within the SPD energy domain, the $y$–distribution serves as a sensitive probe of both the gluon-density normalization and the evolution pattern encoded in different unintegrated PDF parameterizations. The observed energy dependence of the width and normalization provides a consistent phenomenological signature of gluon–gluon fusion dominance and establishes a quantitative baseline for forthcoming SPD measurements of quarkonium production.

It is interesting to note that, at comparable energies, exclusive $J/\psi$ production in proton–proton ultraperipheral collisions yields $d\sigma/dy$ values that are $10^{-5}$–$10^{-4}$ times smaller than those reported here Goloskokov:2025hsk. The large difference originates from the dominance of electromagnetic interactions in ultraperipheral processes. Nevertheless, such measurements provide a valuable complementary perspective on $J/\psi$ production dynamics across different interaction regimes.

The transverse-momentum spectra, $d\sigma/dp_{T}$, presented in Fig. 2, exhibit the expected steep falloff with increasing $p_{T}$. Both gluon-density models reproduce this trend, but the spectral shapes reveal distinct features. At $\sqrt{s}=9~\mathrm{GeV}$, the spectrum is limited to low $p_{T}$, reflecting the restricted phase space near threshold. With $p_{T}$, KL${}^{\prime}2025$ prediction shows a relatively smaller initial amplitude and harder tail near the end. As the energy increases, the spectra broaden from about $p_{T}\!\lesssim\!5~\mathrm{GeV}$ at $\sqrt{s}=9~\mathrm{GeV}$ to nearly $p_{T}\!\sim\!12~\mathrm{GeV}$ at $\sqrt{s}=27~\mathrm{GeV}$, illustrating the gradual transition from a soft, near-threshold domain to a more perturbative regime. Growth in overall normalization is consistent with the increasing gluon luminosity and available phase space. At these energies, LLM${}^{\prime}2024$ prediction relatively yields the initial higher amplitude and harder tail with $p_{T}$ compared to the KL${}^{\prime}2025$ and a crossover region around $2~\mathrm{GeV}<p_{T}<3~\mathrm{GeV}$ is observed, signaling the onset of enhanced partonic activity beyond which KL${}^{\prime}2025$ dominates. These features confirm that the $p_{T}$ spectrum at SPD energies is sensitive to the magnitude and shape of gluon TMDs in the moderate-$x$ domain ($x\sim 10^{-2}\!-\!10^{-1}$).

It is important to mention that LLM’24 density show predicitons beyond the maximum $p_{T}$ range show here but fluctuations in its precision increases rapidly for ${d\sigma}/{dp_{T}}$ values below $10^{-12}~\mathrm{nb/GeV}$.

To further quantify the model differences, normalized $J/\psi$ transverse-momentum spectra are compared in Fig. 3. The data are fitted with empirical functions of the form $A(p_{T}+p_{0})^{n}e^{-Bp_{T}}$ for KL${}^{\prime}2025$ and $A(p_{T}+p_{0})^{n}e^{-Bp_{T}^{2}}$ for LLM${}^{\prime}2024$, yielding excellent agreement across the full $p_{T}$ range.

The KL${}^{\prime}2025$ fit, characterized by a larger parameter $n$, corresponds to a slightly harder rise of the spectrum, whereas the LLM${}^{\prime}2024$ form exhibits a smoother yet faster falloff at high $p_{T}$ due to its Gaussian-like dependence, even though the associated parameter $B$ is smaller. These distinct behaviors directly reflect the different gluon $k_{T}$-broadening mechanisms intrinsic to each TMD parametrization.

Fig. 4 shows the relative contributions of the various $c\bar{c}$ intermediate states to the total $J/\psi$ yield. Both TMD models predict nearly identical channel hierarchies, with CO states—particularly ${}^{3}P_{2}^{[8]}$, ${}^{3}P_{0}^{[8]}$, and ${}^{1}S_{0}^{[8]}$—dominating across the entire SPD energy range. The CS ${}^{3}S_{1}^{[1]}$ contribution remains below the percent level even at $\sqrt{s}=30~\mathrm{GeV}$. This confirms the essential role of CO mechanisms in $J/\psi$ production at intermediate energies.

Finally, the total inclusive cross section as a function of $\sqrt{s}$ is presented in Fig. 5. The cross section rises monotonically with increasing energy and is well described by a shifted power-law form, $(\sigma_{J/\psi}=\alpha(\sqrt{s}-s_{0})^{\beta})$, for $\sqrt{s}\gtrsim 12~\mathrm{GeV}$. The point at $\sqrt{s}=9~\mathrm{GeV}$ lies below the extrapolated curve, suggesting possible near-threshold suppression effects. The LLM${}^{\prime}2024$ parametrization yields a higher normalization and a softer energy dependence ($\beta\simeq 1.8$) compared with KL${}^{\prime}2025$ ($\beta\simeq 2.1$), implying a slower rise of gluon-driven production with energy. Between $12$ and $30~\mathrm{GeV}$, the LLM${}^{\prime}2024$-to-KL${}^{\prime}2025$ scaling ratio decreases from about 2.5 to 2.0, highlighting the diminishing relative contribution from the LLM${}^{\prime}2024$ gluon density as $\sqrt{s}$ increases. These results emphasize the strong sensitivity of $J/\psi$ production to the underlying gluon TMDs in the near-threshold region accessible at SPD/NICA.

Figure 3: Normalized differential cross sections $(1/\sigma)\,d\sigma/dp_{T}$ for $J/\psi$ production in $pp$ collisions at $\sqrt{s}=27~\mathrm{GeV}$ obtained using the KL${}^{\prime}2025$ Kotikov:2025wft and LLM${}^{\prime}2024$ Lipatov:2024xni gluon densities. The solid curves represent empirical fits to the simulated spectra with parameters shown in the legend. The distributions are normalized to emphasize differences in spectral shape between the two TMD sets.

Figure 4: Relative contributions of different intermediate $c\bar{c}$ states to the total $J/\psi$ production cross section in $pp$ collisions as a function of $\sqrt{s}$ for the KL${}^{\prime}2025$ Kotikov:2025wft and LLM${}^{\prime}2024$ Lipatov:2024xni gluon densities. The dominance pattern ${}^{3}P_{2}^{[8]}$$>$${}^{3}P_{0}^{[8]}$$>$${}^{1}S_{0}^{[8]}$$>$${}^{3}P_{1}^{[8]}$$>$${}^{3}S_{1}^{[8]}$$>$${}^{3}S_{1}^{[1]}$ persists across the entire SPD energy range, with minimal model dependence.

Figure 5: Total inclusive $J/\psi$ production cross section in $pp$ collisions as a function of the center-of-mass energy $\sqrt{s}$, obtained with the KL${}^{\prime}2025$ Kotikov:2025wft and LLM${}^{\prime}2024$ Lipatov:2024xni gluon densities. The cross section increases with $\sqrt{s}$, reflecting the growth of gluon luminosity. The shifted power-law fit provides an excellent description for $\sqrt{s}\gtrsim 12~\mathrm{GeV}$; the deviation at $\sqrt{s}=9~\mathrm{GeV}$ likely signals threshold or model-limit effects.

4. Summary— This work presents the first quantitative analysis of inclusive $J/\psi$ production in the SPD/NICA energy regime using modern TMD gluon densities. We have investigated the sensitivity of differential observables to scale variations, intermediate-state contributions, and the collision-energy dependence within both CCFM and KMR-based frameworks. The results establish a robust theoretical baseline for forthcoming SPD experimet and emphasize the need for refined gluon-density extractions in the low to intermediate energy region, where quarkonium production remains a uniquely sensitive probe of nonperturbative gluon dynamics.

Acknowledgements— S.S. expresses sincere gratitude to Prof. A.V. Lipatov for valuable clarifications regarding the use of TMD parton densities in PEGASUS, and to S. Puhan for helpful discussions. This work was done with support from Ministry of Science and Higher Education of Russian Federation under State Assignment № 075-03-2025-662 from 17.01.2025 and by the Committee of Science of the Ministry of Science and Higher Education of the Republic of Kazakhstan (Grant No. BR21881941).