Recent Quarkonia Measurements in Small Systems at RHIC and LHC Energies

Abstract

Heavy-ion research at the Relativistic Heavy Ion Collider (RHIC) during the first decade of data collection, approximately during the years 2000–2010, was primarily focused on the study of Au+Au collisions. The search for evidence of quark-gluon plasma (QGP), a state of matter where quarks and gluons become unbound within a high energy density environment, which was at the forefront of research efforts. However, studies of the azimuthal anisotropy parameter $v_2$ in p/d+Pb collisions from the Large Hadron Collider (LHC) yielded results consistent with the hydrodynamic flow, one of the signatures of quark-gluon plasma formation in heavy-ion collisions. Since the publication of these findings, the field of heavy-ion physics has made subsequent measurements in small system collisions to study cold nuclear matter effects as well as look for additional evidence of hot nuclear matter effects. Quarkonia, a bound state of a $c\overline{c}$ or $b\overline{b}$ pair, has often been used to probe a wide range of nuclear effects in both large and small collision systems. Here we will review recent quarkonia measurements in small system collisions at RHIC and LHC energies and summarize the experimental conclusions.

1. Introduction

Right around the time of the development of quantum chromodynamics in the 1970s, Brookhaven National Laboratory (BNL) submitted to Physical Review Letters (PRL) evidence for the discovery of a new particle, a meson made up of a charm and an anti-charm quark [1]. The next day, Stanford Linear Accelerator submitted to PRL evidence for their discovery of the very same particle [2]. Both experimental results were published on 2 December 1974, in the same issue of PRL. Two years later, Professors Samuel C.C. Ting of BNL and Burton Richter of SLAC received the Nobel Prize for discovering this new particle, now known as the $J/\psi$ meson (see Figure 1).

Meanwhile, across the Atlantic Ocean in the offices of the European Organization for Nuclear Research (CERN), theoretical physicist Rolf Hagedorn had been working for years on the statistical thermodynamics of strong interactions at high energies. In 1965, he discovered there existed a maximum temperature for the strong interaction, which he determined to be $T_H\sim$ 158 MeV [3]. However, after the introduction of asymptotic freedom [4] by David Gross and Frank Wilczek, Hagedorn came to understand that the Hagedorn temperature was less of a limiting temperature than a reflection of the phase transition from a hadron gas, confined by the strong interaction, into a quark-gluon plasma, deconfined by asymptotic freedom.

By 1980, the International Conference on Ultrarelativistic Nucleus–Nucleus Collisions (Quark Matter) series had begun. And by 1984, a formal proposal had been put forth to build the Relativistic Heavy-Ion Collider (RHIC) at BNL [5]. Combining the discoveries of the $J/\psi$ meson and asymptotic freedom, Tetsuo Matsui of the Massachusetts Institute of Technology and Helmut Satz of Bielefeld University published a paper in the 1980s predicting an “unambiguous signature” for the formation of quark-gluon plasma in heavy-ion collisions. This signature was indeed the suppression of $J/\psi$ production in nuclear–nuclear (A+A) collisions with respect to $J/\psi$ production in proton-proton (p+p) collisions as a result of color screening in a deconfined medium.

CERN was already running ion beams O and S with the NA38 experiment at the Super Proton Synchrotron (SPS) in 1989. A few years later, the NA50 experiment started collecting PbPb interaction data. Designed to measure muons from vector meson decays, NA50 released the $J/\psi$ cross section measurements (Scaled by the factor $AB$, the product of the projectile A and target B mass numbers) from NA38 p+A collisions and NA51 AB collisions to compare with the new $J/\psi$ cross section measurements (Scaled by the factor $AB$, the product of the projectile A and target B mass numbers) from Pb+Pb collisions [6]. In Figure 2, the $J/\psi$ measurements taken in p+A collisions fall along a similar trendline. But the $J/\psi$ suppression observed in Pb+Pb collisions was ∼5 standard deviations below the trend seen in p+A collisions, with the behavior described as anomalous suppression.

RHIC at BNL completed construction in the late 1990s and started recording its first data in 2000, where the “main goal of the RHIC heavy ion program is the discovery of the novel ultra–hot, high–density state of matter predicted by the fundamental theory of strong interactions and created in collisions of heavy nuclei, the Quark–Gluon Plasma (QGP) [7]”. In 2005, all four experimental programs at RHIC at the time—PHENIX, STAR, BRAHMS, and PHOBOS—announced the formation of dense partonic matter in Au+Au collisions consistent with quark-gluon plasma [8,9,10,11]. Results from Pb+Pb data were available from the LHC as soon as 2011 [12,13], and the collective behavior of the highly dense system was modeled using viscous fluid dynamics in a 2012 paper by theorists Charles Gale, Sangyong Jeon, Björn Schenke, Prithwish Tribedy, and Raju Venugopalan [14]. They showed that hydrodynamics successfully describes the anisotropic flow coefficients measured by ATLAS and PHENIX (see Figure 3), as well as for ALICE and STAR.

However, earlier that year, in 2012, theorist Piotr Bożek of the Rzeszów University in Poland published a paper entitled, “Collective flow in p-Pb and d-Pb collisions at TeV energies” [15]. This publication was essentially the first time in the era of heavy-ion physics that the idea of quark-gluon plasma formation in small collision systems had ever been suggested. Although the wider heavy ion community had assumed QGP was not formed in small system collisions, in 1983, Rolf Hagedorn believed it was already being created at CERN’s Internal Storage Ring (ISR) energies [16]. In Hagedorn’s concluding statements, he asks the question, “Do we see the phase transition hadron → quark-gluon plasma (predicted by so many models) at $p\overline{p}$ collider energies? Yes, we even see it already at ISR energies [16]”.

2. Hot versus Cold Nuclear Matter Effects

For more clarity, when discussing the quarkonia results, we will define the commonly used terms and what is meant by hot and cold nuclear matter effects. Cold nuclear matter effects occur independently of quark–gluon plasma formation in the system [17]. These effects include initial-state energy loss [18], gluon shadowing and anti-shadowing [19,20], nuclear absorption (or nuclear break-up) [21,22], the co-mover particle interaction model [23,24], and the Cronin effect [25]. On the other hand, hot nuclear matter effects generally refer to quark-gluon plasma formation as well as regeneration (or coalescence) [26]. Cold nuclear matter effects are expected to occur alongside hot nuclear matter effects in A+A collisions.

Here we will briefly describe each effect, but for more information, please see the cited references. Initial-state energy loss refers to the projectile gluon undergoing multiple scattering while passing through the target before $J/\psi$ production. The co-mover interaction model (CIM), first introduced in the 1990s [27], involves no phase transition and instead describes the $J/\psi$ suppression as a result of break-up from co-moving particles created during the collision. Gluon shadowing and anti-shadowing result from modifications to the nucleus’s gluon nuclear parton distribution functions (nPDFs). And lastly, nuclear absorption refers to the break-up of the bound $J/\psi$ (or precursor state) in collisions with other target nucleons passing through the $J/\psi$ production point.

3. Quarkonia at the Relativistic Heavy-Ion Collider

We review recent quarkonia results from the PHENIX and STAR Collaborations at the Relativistic Heavy Ion Collider. Currently, STAR is still recording data, with beam time scheduled through the end of 2025. PHENIX, however, ended its data-taking run in 2016 but is still analyzing what is considered its best data sets on record (Runs 14, 15, and 16).

3.1. PHENIX Collaboration Results

The nuclear modification of $J/\psi$ in Au+Au collisions was an early measurement published by the PHENIX Collaboration [28], which started recording data in the year 2000. As the years progressed, upgrades were installed, and the recorded luminosity increased. Shown on the left of Figure 4 [29] is the PHENIX $J/\psi$ nuclear modification from Run 4 (blue data points) [30] as a function of multiplicity density compared with the same measurement made with the ALICE detector (red data points) [31]. Consistent with the predictions by Matsui and Satz and similar to the NA50 Collaboration findings, the $J/\psi$ suppression observed in Au+Au collisions by PHENIX is extremely strong. The results at the lower RHIC energies of 200 GeV per nucleon-nucleon collision are consistent with the final state effect of quark–gluon plasma formation. However, we see distinctly different behavior at the LHC, where the nuclear modification does not reach the level of suppression seen at the RHIC. This striking difference can be explained by regeneration, which can modify $J/\psi$ yields at the higher LHC energies.

On the right of Figure 4, a theoretical prediction of regeneration (or coalescence) at LHC energies by Du & Rapp [32] is shown. The blue dashed line represents the contribution to the $J/\psi$ nuclear modification factor due to regeneration, while the orange dashed line represents the total contribution from primordial $J/\psi$. Neglecting the contributions from recombination, we can see the predicted $J/\psi$ nuclear modification factor at LHC energies is likely more suppressed than what we observe at RHIC energies. However, due to the regeneration effects at the LHC, which are not present at RHIC energies, the total $J/\psi$ nuclear modification factor is less suppressed at ALICE than at PHENIX. Nonetheless, the results from ALICE are also consistent with the final state effect quark-gluon plasma formation.

The $J/\psi$ nuclear modification factor (filled data points) in d+Au collisions is shown in Figure 5, compared with muons from open heavy flavor decays (open data points) [33]. All data are from the 0–20% centrality class and the red (blue) data points represent forward (backward) rapidity. At forward rapidity, the $J/\psi$ suppression is quite similar to the open charm suppression. At RHIC energies, the open heavy flavor muons are predominately from D-meson decays. Suppression would be consistent with shadowing and/or initial state parton energy loss, similarly modifying the $J/\psi$ and open heavy flavor production.

At backward rapidity, the $J/\psi$ is suppressed relative to the open charm. It is expected that open charm should be enhanced by anti-shadowing, as the data indicate, with the nuclear modification factor well above unity. However, the suppression of the $J/\psi$ is consistent with nuclear absorption due to collisions with other nucleons inside the Au target. Previous publications have also noted the possibility of nuclear absorption [34,35,36,37,38,39,40], although it is possible that there could be a contribution from comovers as well.

In Figure 6, the inclusive $J/\psi$ nuclear modification factor is shown as a function of $p_T$ in p+Au collisions at forward (left) and backward rapidity for 0–100% centrality [41]. Suppression is observed in the forward rapidity region for approximately $p_T$ < 3 GeV/c. Above approximately $p_T$ > 3 GeV/c, the data slowly increases until the modification is above unity. A similar trend is seen in the data at backward rapidity, with suppression at lower $p_T$ and rising to around or above unity. The suppression is not as pronounced and is only observed for approximately $p_T$ < 2 GeV/c.

The measurements are compared with gluon nPDF predictions provided by Shao et al. [42,43,44,45], Du & Rapp [26,46,47,48], and R. Vogt [49,50] (see [41] for model descriptions). The striped magenta curves are Vogt EPPS16 Leading-Order LO predictions. The solid blue (gray lattice) curves are Shao re-weighted EPPS16 (nCTEQ15) predictions, and lastly, the red lattice curves are Du & Rapp transport model predictions. At forward rapidity, all models describe the data well. The suppression seen in the transport model is primarily due to gluon shadowing. The data is slightly more suppressed than the transport model predicts, possibly due to the older EPS09 nPDFs. Both the re-weighted EPPS16 and re-weighted nCTEQ15 predictions (which include LHC data) describe modification well and have significantly improved uncertainty over the EPPS16 NLO predictions. At backward rapidity, anti-shadowing effects alone from EPPS16 and nCTEQ15 predictions cannot describe the suppression seen in the data at low $p_T$. The $J/\psi$ modification at backward rapidity is best described by the transport model, which also includes a nuclear absorption estimate.

In Figure 7, a comparison of inclusive $J/\psi$ and $\psi$(2S) nuclear modification is shown at both RHIC and LHC energies [51]. The blue (gold) data points are PHENIX measurements in p+Au (d+Au) collisions. The red (gray) data points are LHCb (ALICE) measurements in p+Pb collisions. The open (filled) data points represent the $J/\psi$ ($\psi$(2S)) nuclear modification measurements. Cold nuclear matter effects are expected to be similar for the two charmonium states, so we would expect nuclear absorption and gluon shadowing/anti-shadowing to be similar as long as the collision system, collision energy, and rapidity range are the same for the $J/\psi$ and $\psi$(2S) measurements. We can see similar suppression between the two states at forward rapidity, consistent with gluon shadowing as the dominant contribution. At backward rapidity, however, the $\psi$(2S) is much more suppressed than the $J/\psi$, indicating that the $\psi$(2S) is more susceptible to final state effects than the $J/\psi$. Those final state effects are generally considered either co-moving particle interactions or quark–gluon plasma formation.

We note that at the time of writing, preliminary results for $J/\psi$ elliptic flow ($v_2$) in Au+Au collisions and multiplicity-dependent $J/\psi$ and $\psi$(2S) production in p+p collisions are available, but the studies have not been formally published, and so we omit these results from this review.

3.2. STAR Collaboration Results

The STAR Collaboration has recently measured the $J/\psi$ nuclear modification factor through the dimuon decay channel as a function of $N_{part}$ in Au+Au collisions [52]. We include these results (shown in Figure 8) as a comparison to the PHENIX $J/\psi$ results in Au+Au collisions shown in Figure 4. In the left panel of Figure 8, the $J/\psi$ nuclear modification results are essentially $p_T$ integrated ($p_T$ > 0.15 GeV/c for STAR, and $p_T$ > 0 GeV/c for ALICE). The modification in the right-hand panel is shown at $p_T$ > 5 Gev/c for STAR measurements and $p_T$ > 6.5 GeV/c for ALICE measurements. Charm quark coalescence (or recombination) is expected to occur at low $p_T$, likely below $p_T$ < 6.5 GeV/c.

In the left-hand plot, the STAR $J/\psi$ suppression is stronger than the suppression observed at ALICE, consistent with the previously shown PHENIX results. Note that the STAR measurements were taken at mid-rapidity, while the PHENIX measurements were taken at forward rapidity. However, in the right-hand plot for the higher $p_T$ range, the STAR data points fall slightly above the ALICE data points. The comparison at higher $p_T$ is expected to exclude cold nuclear matter effects. Since the suppression is slightly stronger at LHC energies, this could indicate that the temperature of the nuclear medium created in Pb+Pb collisions at LHC energies ($\sqrt{s_{NN}}=2.76$ TeV) is higher than that created in Au+Au collisions at RHIC energies ($\sqrt{s_{NN}}=0.2$ TeV).

The data are compared with three different theoretical models (see [52] for model descriptions). In both the left and right plots, the magenta dashed line and the gray lattice curve represent transport model (TM) predictions from Tsinghua University [53] and Texas A&M University [46], respectively. In the left-hand plot only, there is an additional comparison to the statistical hadronization model (SHM). Both transport model predictions describe the data well at low $p_T$ in the left-hand figure. In contrast, the statistical hadronization model describes the data well up until the peripheral events (∼$N_{part}<100$). At higher $p_T$ in the right-hand figure, the transport model predictions appear to either over-predict or under-predict the STAR $J/\psi$ suppression.

In Figure 9, the inclusive $J/\psi$ polarization parameters ${\lambda }_{\theta }$, ${\lambda }_{\varphi }$, and ${\lambda }_{\theta \varphi }$ are shown in p+p collisions at mid-rapidity as a function of $p_T$ [54]. The open (filled) data points represent the $J/\psi$$\to e^{+}{e}^{-}$ ($J/\psi$$\to {\mu }^{+}{\mu }^{-}$) decay channel. The left column shows the polarization parameters measured in the Helicity reference frame (HX), and the right column shows the parameters measured in the Collins-Soper reference frame (CS). We can see that nearly all parameters in both frames of reference are consistent with zero within uncertainties, except for the largest $p_T$ measurement for the ${\lambda }_{\theta }$ parameter in the Collins-Soper frame. However, overall, no meaningful transverse (+1) or longitudinal (−1) $J/\psi$ polarization can be determined from the measurements.

The data are compared with charmonium production models (see [54] for model descriptions). Although the theoretical model predictions are for prompt $J/\psi$ production, the contribution of non-prompt $J/\psi$ from B-decays in the inclusive sample is expected to be small. The pink (cyan) striped curves denote the next-to-leading-order non-relativistic quantum chromodynamics NLO NRQCD1 [55] (NLO NRQCD2 [56]) prediction models for the rapidity range $\left|y\right|<1.0$, while the red (blue) denote the same models for the rapidity range $\left|y\right|<0.5$. The turquoise (green) striped curves represent the color glass condensate + NRQCD [57] for a rapidity range $\left|y\right|<1.0$ ($\left|y\right|<0.5$). And the orange curve denotes the improved color evaporation model [58]. Although none of the models can be decisively ruled out, the color glass condensate + NRQCD best describes the data overall.

In Figure 10, the $J/\psi$ production in p+p collisions is shown as a function of multiplicity [59]. Along the x-axis, the charged particle multiplicity per pseudo-rapidity $d{N}_{ch}^{MB}/d\eta$ is scaled by the average charged particle multiplicity per pseudo-rapidity $\langle d{N}_{ch}^{MB}/d\eta \rangle$, and along the y-axis, the $J/\psi$ yield $N_{J/\psi }$ is scaled by the average $J/\psi$ yield $\langle N_{J/\psi }\rangle$. The open data points in the left panel represent ALICE measurements at $p_T$ > 0. The purple data points represent measurements at $p_T$ > 0 GeV/c, while the blue (red) data points represent measurements at $p_T$ > 1.5 ($p_T$ > 4) GeV/c. The relative $J/\psi$ yield increases as the relative charged particle multiplicity increases for all $p_T$ ranges. The measurements by ALICE ($\sqrt{s}=7$ TeV) and STAR ($\sqrt{s}=0.2$ TeV) are surprisingly similar, despite an order of magnitude difference in collision energy. This similarity suggests that the underlying mechanism for $J/\psi$ production is independent of collision energy.

The data are compared with Monte Carlo generator calculations (see [54] for more details and MC references). The middle panel of Figure 10 shows three curves generated by PYTHIA8 [60] that correspond to each $p_T$ range of the data. Multi-partonic interactions, where more than one interaction per p+p collision occurs, were included in the PYTHIA calculation and can describe the increasing trend seen in the data as a function of charged particle multiplicity. In the far right panel, the data are compared with EPOS3 [61] calculations at higher $p_T$ and the percolation model [62] for $p_T>0$ GeV/c. Note that the EPOS3 curves correspond to open charm production and slightly different $p_T$ ranges than the data. The EPOS3 calculations also include multi-partonic interactions and agree well with the data. The percolation model, although developed for high-density p+p collisions at LHC energies, describes the STAR data well at $p_T>0$ GeV/c.

Figure 11 shows the $J/\psi$ nuclear modification factor in p+Au, d+Au, and Au+Au collision systems as a function of $p_T$ at mid-rapidity [63]. The red (blue) data points represent STAR measurements in p+Au (Au+Au) collision systems, while the open data points represent PHENIX measurements in d+Au collisions. The $J/\psi$ suppression as a function of $p_T$ in Au+Au collisions is flat, without any noticeable $p_T$ dependence. The small system measurements in p+Au and d+Au collisions, however, do show a strong $p_T$ dependence. It is also worth mentioning that the STAR and PHENIX measurements are quite similar, particularly the suppression at $p_T$ < 4 GeV/c, despite the different projectile in each collision system. Note that all measurements were made at mid-rapidity. The suppression seen above 3 GeV/c is expected to be outside the threshold for cold nuclear matter effects and could indicate quark-gluon plasma formation modifying the $J/\psi$ production in small system p+Au collisions.

The data are compared with several different theoretical models (see [63] for model descriptions), including the improved color evaporation model (ICEM [64]) and the color glass condensate (CGC [65]). The beige (navy) curves represent the ICEM+NLO EPS09 (CGC+ICEM) predictions. The light green (olive green) curves represent the TAMU [47,48] (Lansberg [42,43,44,45]) predictions, and the dashed purple (solid magenta) lines represent the comover interaction [24] (energy loss + broadening [66]) models. The comover interaction model shows little dependence on the $J/\psi$ transverse momentum, predicting a similar suppression at low and high $p_T$, as observed in Au+Au collision system data. All theoretical models reasonably describe the data and reproduce the suppression at low transverse momentum.

We note that at the time of writing, preliminary results for $\mathsf{{\rm Y}}$(1S) and $\mathsf{{\rm Y}}$(2S) nuclear modification in Au+Au collisions are available [67], but the study has not been formally published and so we omit these results from this review.

4. Quarkonia at the Large Hadron Collider

We review recent quarkonia results from the ALICE, ATLAS, CMS, and LHb Collaborations at the Large Hadron Collider. All four Collaborations have been actively recording data since the year 2010. In July 2022, the collaborations began Run 3 of data taking. Due to rising energy costs, a 20% reduction in LHC operation time has been announced for the 2023 schedule.

4.1. ALICE Collaboration Results

The $J/\psi$ nuclear modification is shown in Figure 12 as a function of $p_T$ at forward rapidity (left) and backward rapidity (right) in p+Pb collisions for the 2–10% centrality class [68]. The red (blue) data points correspond to $\sqrt{s_{NN}}=8.16$ (5.02) TeV. The modification appears independent of collision energy, as both sets of measurements are consistent. At forward rapidity, suppression is observed for approximately $p_T$ < 6 GeV/c, where cold nuclear matter effects are expected to dominate. The $J/\psi$ nuclear modification dependence on $p_T$ is similar to the results seen in p+Au collisions at PHENIX in Figure 6 and at STAR in Figure 11. However, the suppression at RHIC energies is observed for approximately $p_T$ < 3 GeV/c. At backward rapidity, the $p_T$ dependence in p+Pb collisions looks different from the $p_T$ dependence at forward rapidity. Minimal suppression is observed, but only for around $p_T$ < 1 GeV/c. The nuclear modification measurements for around $p_T$ > 1 GeV/c are more or less above unity. This lack of suppression at low $p_T$ contrasts with the PHENIX results at backward rapidity shown in Figure 6, where the dependence is similar to what is observed at forward rapidity. This contrast suggests the dominant source of $J/\psi$ modification at backward rapidity differs at RHIC and LHC energies.

The data are compared with several different theoretical models (see [68] for model descriptions). The striped red curve corresponds to the Vogt et al. EPS09s NLO + color evaporation model [21]. The solid green curve corresponds to the Arleo and Peignè energy loss model [69], and the striped purple curve corresponds to the Du & Rapp transport model [47]. At forward rapidity, the data is best described by the Color Evaporation and transport models. At backward rapidity, the data shows strong anti-shadowing for the most central collisions. Each model predicts some suppression, with the energy loss and transport models predicting the most significant degree of suppression at low $p_T$. However, the modification is squarely above unity for approximately $p_T$ > 2 GeV/c. There does not appear to be any $J/\psi$ modification due to final state effects at backward rapidity in central p+Pb collisions at ALICE.

In Figure 13, the $\mathsf{{\rm Y}}$(2S) to $\mathsf{{\rm Y}}$(1S) ratio of yields (left) and $R_{AA}$ (right) in Pb+Pb collisions are shown at forward rapidity ($2.5<y<4.0$) as a function of $\langle N_{part}\rangle$ [70]. The ratio at higher $\langle N_{part}\rangle$ and more central collisions is more suppressed than at the lower $\langle N_{part}\rangle$ and more peripheral collisions, as would be expected due to final state effects.

The data are compared with several different theoretical models (see [70] for model descriptions). In the left-hand figure, the ratio of $\mathsf{{\rm Y}}$(2S) to $\mathsf{{\rm Y}}$(1S) yields is compared with the statistical hadronization model [71], which overestimates the suppression of the $\psi$(2S) relative to the $J/\psi$. In the right-hand figure, the ratio of the $\mathsf{{\rm Y}}$(2S) to $\mathsf{{\rm Y}}$(1S) nuclear modification is compared with hydrodynamics calculations [72], transport model predictions with and without regeneration effects [73], calculations based on coupled Boltzmann equations [74], and the comover interaction model [23]. The data is not well described by the comover interaction model. The hydrodynamic calculations and the transport model with coalescence effects are most consistent with the measured data, indicating the final state effects are more consistent with quark-gluon plasma formation.

Figure 14 shows the $J/\psi$ elliptic flow $v_2$ as a function of $p_T$ in Pb+Pb collisions for the 20–40% centrality class [75]. The red (blue) data points represent the $J/\psi$$\to e^{+}{e}^{-}$ ($J/\psi$$\to {\mu }^{+}{\mu }^{-}$) decay channels. The elliptic flow measurements are unambiguously positive, with the $v_2$ coefficient for the 4–6 GeV/c $p_T$ range carrying a 6.6$\sigma$ statistical significance. These results clearly show that $J/\psi$ mesons participate in collective behavior.

The data are compared with several different transport model calculations (see [75] for model descriptions). The blue (red) lattice curve represents the Du & Rapp transport model [47] at mid (forward) rapidity, and the red dashed line is the contribution from primordial $J/\psi$. The orange open curve (dashed line) represents the K. Zhou et al. transport model [76] at forward rapidity with (without) non-collective flow, while the orange dotted line corresponds to the primordial $J/\psi$ contribution. The transport model calculations describe the $v_2$ measurements well at mid-rapidity but do not entirely capture the behavior for forward rapidity data above approximately $p_T$ > 4 GeV/c. The K. Zhou et al. transport model predicts a larger contribution of elliptic flow at high $p_T$ but again does not entirely capture the collective behavior. It has been suggested that magnetic field effects could potentially be missing from the transport model calculations [75].

We note that at the time of writing, preliminary results for the $\psi$(2S) to $J/\psi$ ratio in Pb+Pb collisions [77], $J/\psi$ production at mid-rapidity in p+Pb collisions [78], and the $J/\psi$ polarization in Pb+Pb collisions [79] are available, but the studies have not been formally published, and so we omit these results from this review.

4.2. ATLAS Collaboration Results

The elliptic flow $v_2$ of muons from p+p collisions at $\sqrt{s}=$ 13 TeV is shown in Figure 15, where charm decays are represented by the blue data points and bottom decays by the red data points [80]. In the left-hand figure, $v_2$ is plotted as a function of charged track multiplicity $N_{ch}^{rec}$ and as a function of $p_T$ in the right-hand figure. The elliptic flow of muons from the bottom decays independently of the charge track multiplicity and is consistent with zero within uncertainties. Similarly, the muon $v_2$ from bottom quark decays shows no dependence on $p_T$ and is again consistent with zero within uncertainties. In contrast, the muons from the lighter charm quark decay with non-zero elliptic flow as a function of charge track multiplicity and as a function of $p_T$. These results indicate collective behavior in small system collisions and that while the lighter charm quarks participate, the heavier bottom quarks do not.

The $\psi$(2S) to $J/\psi$ ratio in p+Pb collisions is shown in Figure 16 in the mid-rapidity region $-2<y^{*}<1.5$ as a function of centrality [81]. The measurements were taken over the transverse momentum range $8<p_T<40$ GeV/c. It is expected that nuclear effects related to the initial state will cancel by taking the ratio between two charmonium states. These effects could include initial state energy loss or gluon shadowing, which depend on the target nucleus, collision energy, and rapidity and, therefore, would be the same for both states. The ratio’s final state effects remain, including the comover interaction, quark-gluon plasma formation, or potentially nuclear absorption effects (although nuclear absorption is not believed to make an impact at LHC energies due to the large Lorentz factor). The ratio of $\psi$(2S) to $J/\psi$ is right around unity for peripheral collisions, and no clear difference in modification between the two states is observed. However, the $\psi$(2S) is suppressed with respect to the $J/\psi$ for the more central 0–5%, 5–20%, and 20–40% centrality classes. Although coalescence is not expected to significantly contribute to charmonium production in p+Pb collisions, the $p_T$ range of the measurements at $8<p_T<40$ GeV/c likely falls beyond the region where coalescence effects are believed to occur. The uncertainties prevent any firm conclusions, but the data suggest the $\psi$(2S) is more suppressed than the $J/\psi$ in central p+Pb collisions due to final state effects. These ATLAS results are consistent with the PHENIX and LHC data shown in Figure 7, where stronger suppression is observed at backward rapidity for $\psi$(2S) with respect to $J/\psi$ due to final state effects in p+A collisions.

We note that, at the time of writing, preliminary results for centrality-dependent jet quenching in p+Pb collisions are available [82], but the studies have not been formally published, and so we omit these results from this review.

4.3. CMS Collaboration Results

In Figure 17, the elliptic flow $v_2$ is shown as a function of $p_T$ for the $D^0$ and $J/\psi$ mesons in p+Pb collisions [83]. The $D^0$ $v_2$ measurements (red data points) and $J/\psi$ $v_2$ measurements (blue data points) are compared with the $\mathsf{\Lambda }$ baryon (open circles) and $K_S^0$ meson (open squares) measurements. These measurements are for prompt production, where the $D^0$ and $J/\psi$ mesons were produced in the primary interaction and would therefore experience the full effect of the nuclear medium if one were created. The $K_S^0$ and $\mathsf{\Lambda }$ hadrons are composed of light and strange-flavored quarks and show a larger elliptic flow than the $J/\psi$ and $D^0$ mesons that contain charm. A clearly positive $v_2$ is observed for both the $D^0$ and $J/\psi$ mesons. Consistent with the results from ATLAS shown in Figure 15, CMS also finds evidence of collective flow in small system collisions. These CMS results support the ATLAS conclusion that charm quarks participate in collective behavior in small system collisions.

Figure 18 shows the nuclear modification of $\mathsf{{\rm Y}}\left(1S\right)$, $\mathsf{{\rm Y}}\left(2S\right)$, and $\mathsf{{\rm Y}}\left(3S\right)$ in p+Pb collisions as a function of rapidity [84]. The left-hand plot shows the bottomonium measurements for transverse momentum $p_T$ < 6 GeV/c, and the right-hand plot shows the measurements for the range $6<p_T<30$ GeV/c. The forward rapidity ($0<y<1.93$) measurements for the three states do not significantly change between the two $p_T$ ranges, which likely indicates cold nuclear matter effects are the dominant source of modification. At backward rapidity ($-1.93<y<0$), the suppression of the $\mathsf{{\rm Y}}\left(3S\right)$ is very strong compared with the $\mathsf{{\rm Y}}\left(1S\right)$ at low $p_T$. The suppression between the two states is pronounced, although not as dramatic in the higher $p_T$ range. It is interesting to note at low $p_T$ that the $\mathsf{{\rm Y}}\left(1S\right)$ is less suppressed at backward than forward rapidity, the $\mathsf{{\rm Y}}\left(2S\right)$ is approximately equally suppressed, and the $\mathsf{{\rm Y}}\left(3S\right)$ is noticeably more suppressed. This sequential suppression of the three bottomonium states supports final state effects in p+Pb collisions.

4.4. LHCb Collaboration Results

In Figure 19, the prompt $J/\psi$ nuclear modification factor (left) is compared with the nonprompt $J/\psi$ nuclear modification (right) as a function of rapidity in p+Pb collisions [85]. For the prompt $J/\psi$, the nuclear modification is flat and does not show a strong dependence on rapidity, unlike the PHENIX inclusive $J/\psi$ results from p+Au collisions shown in Figure 6. The suppression is stronger at forward rapidity ($1.5<y<4.0$) than at backward rapidity ($-5.0<y<2.5$) for prompt production, which is also predicted by the models. The nuclear modification shows slight suppression at backward rapidity, with all measurements being less than unity. This suppression contrasts with the ALICE measurements shown in Figure 12, where the modification is generally larger than unity at backward rapidity. This difference could possibly be due to the $J/\psi$ sample measured: ALICE measured inclusive $J/\psi$, while LHCb measured prompt $J/\psi$.

The data are compared with charmonium modification models, including gluon shadowing and energy loss (see [85] for model descriptions). The beige (blue) curve represents the EPS09 LO [86] (EPS09 NLO [87]) nuclear parton distribution functions, while the dashed orange (solid violet) curve represents the nDSg LO parameterization [86] (fully coherent energy loss FCEL model [88,89]). At forward rapidity, suppression of the $J/\psi$ nuclear modification factor is observed, consistent with RHIC and LHC data shown in Figure 7. The prompt contribution shows more substantial suppression than the nonprompt contribution. In both cases, the suppression is well described by the EPS09 LO gluon shadowing predictions. The nDSg LO model does not predict strong enough suppression for the prompt $J/\psi$ contribution; however, it describes the nonprompt $J/\psi$ measurements well. At backward rapidity, the EPS09 models predict stronger anti-shadowing effects than what is seen in the data. The nDSg LO and FCEL models best describe the hint of suppression at backward rapidity, suggesting the possibility of a nuclear effect beyond gluon anti-shadowing. The uncertainties are more considerable in the backward rapidity measurements for the nonprompt contribution, and both nDSg LO and EPS09 LO describe the data reasonably well.

The ratio of ${\chi }_{c1}$(3872) to $\psi$(2S) from prompt (red squares) and b decays (black squares) as a function of charged tracks in p+p collisions [90] is shown in Figure 20. The prompt contribution generally follows a decreasing trend with increasing multiplicity, suggesting that the ${\chi }_{c1}$(3872) is more suppressed with respect to the $\psi$(2S). The slope of the decreasing measurements is five standard deviations away from a zero slope. The contribution from b decays has the opposite behavior, following an increasing trend with larger multiplicity.

The data are compared with different post-diction comover interaction models [24,91,92] (see [90] for model descriptions). The gray (blue) curves represent the molecular interpretation of ${\chi }_{c1}$(3872) based on coalescence (geometry), while the red curve represents the compact tetraquark interpretation. The compact tetraquark interpretation provides the best description of the data. The molecular interpretations have a much stronger dependence on multiplicity and either immediately fall or rise with an increasing number of charged tracks.

In Figure 21, the $\mathsf{{\rm Y}}$(1S) (left) and $\mathsf{{\rm Y}}$(2S) (right) nuclear modifications are shown as a function of rapidity in p+Pb collisions [93]. Both bottomonium states show suppression at forward rapidity. Although there is some scatter in the $\mathsf{{\rm Y}}$(2S) measurements, the two states appear to experience a similar degree of medication. At backward rapidity, the $\mathsf{{\rm Y}}$(2S) measurements are all below unity, showing a slightly stronger suppression with respect to the $\mathsf{{\rm Y}}$(1S). For comparison, the binding energy of the $J/\psi$ ($\psi$(2S)) is 0.64 (0.05) GeV, while the binding energy of the $\mathsf{{\rm Y}}$(1S) ($\mathsf{{\rm Y}}$(2S)) is 1.1 (0.53) GeV [94], such that the binding energy can be ordered from lowest to highest as $\psi$(2S) < $\mathsf{{\rm Y}}$(2S) < $J/\psi$ < $\mathsf{{\rm Y}}$(1S). Based on binding energy alone, $\mathsf{{\rm Y}}$(1S) is not expected to show a similar degree of suppression as the $J/\psi$. However, there could be an additional influence related to bottomonium formation time compared with charmonium formation time.

The data are compared with different bottomonium modification models (see [93] for model descriptions). The beige (blue) curves represent the EPPS16 [19,45] (nCTEQ15 [20,43]) nuclear parton distribution functions, while the open green (open magenta) curves represent EPS09LO (nCTEQ15) plus the comovers interaction model [23]. At forward rapidity, there is no significant difference between the four models, and all described both the $\mathsf{{\rm Y}}$(1S) and $\mathsf{{\rm Y}}$(2S) modifications well. These results are again compatible with gluon shadowing as the dominant source of suppression at forward rapidity, consistent with other quarkonia results from RHIC and LHC energies. The results are slightly less clear at backward rapidity with the $\mathsf{{\rm Y}}$(1S) state due to the scatter, although each model partially describes the data. The $\mathsf{{\rm Y}}$(2S) results at backward rapidity are more suppressed than what is predicted by anti-shadowing effects alone. Including the comover interaction model provides the best description of the data.

We note that at the time of writing, preliminary results for the $J/\psi$ to $D^0$ ratio using the LHCb SMOG (System for Measuring Overlap with Gas) detector is available [95], but the study has not been formally published and so we omit these results from this review.

5. Discussion

Strong $J/\psi$ suppression is undoubtedly observed in A+A collisions at RHIC (Figure 4) and SPS energies (Figure 2), as Matsui and Satz predicted back in 1986. However, this same degree of suppression is not observed at LHC energies (Figure 4). The (likely) higher energy density of the nuclear medium produced at the Large Hadron Collider leads to a larger fraction of charm quark thermalization and elliptic flow (Figure 14). Combined with a higher $c\overline{c}$ pair production cross section, charm quark coalescence is more significant at LHC versus SPS and RHIC energies. The measurement of the $\mathsf{{\rm Y}}$(2S) to $\mathsf{{\rm Y}}$(1S) ratio in Pb+Pb collisions (Figure 13) shows $\mathsf{{\rm Y}}$(2S) suppression is best described by thermodynamic calculations and Transport Model predictions, supporting the presence of a hot and dense medium formed in ultra-relativistic heavy ion collisions.

In small system collisions, different nuclear effects dominate, as shown by the comparison between $J/\psi$ nuclear modification in Au+Au versus d+Au collisions (Figure 11). At forward rapidity, all quarkonia measurements at RHIC and the LHC show suppression consistent with cold nuclear matter effects. These measurements include $J/\psi$ at PHENIX (Figure 6), charmonia at PHENIX, ALICE, and LHCb (Figure 7), $J/\psi$ at ALICE (Figure 12), and bottomonium at both CMS (Figure 18) and LHCb (Figure 21). ALICE and PHENIX see similar $J/\psi$ modification as a function of $p_T$, with the strongest suppression at low $p_T$. This similarity suggests the source of suppression could be the same at RHIC and LHC energies. None of the models presented, including gluon shadowing predictions, the Transport Model, the Color Evaporation Model, or the comover interaction model, can be definitively ruled out. However, the models include gluon shadowing as the primary contribution, and can include additional nuclear effects.

At backward rapidity, the RHIC and LHC experiments probe the anti-shadowing region of Bjorken-x, where in the absence of other nuclear effects, an enhancement is expected in quarkonia nuclear modification measurements. $J/\psi$ results from LHCb (Figure 19), ALICE (Figure 12), and PHENIX (Figure 6) all show different modifications. The LHCb measurements for prompt $J/\psi$ nuclear modification show slight suppression at backward rapidity. The ALICE inclusive $J/\psi$ measurements show essentially no suppression with unambiguous enhancement, while the PHENIX inclusive $J/\psi$ results show the strongest suppression of the three experiments. Additionally, the ALICE and PHENIX $J/\psi$ modifications as a function of $p_T$ look significantly different, suggesting that nuclear effects on $J/\psi$ production at backward rapidity are different at RHIC and LHC energies. The PHENIX open heavy flavor results (Figure 5) support the conclusion that nuclear absorption effects are significant at RHIC energies. It is possible that the LHCb prompt $J/\psi$ results reflect either hot nuclear matter effects or break-up from comover interactions that are not present in the inclusive $J/\psi$ measurements from ALICE. Transport Model predictions, including hot nuclear matter and nuclear absorption effects, describe the PHENIX data well at backward rapidity. From this comparison of results, it appears the ALICE measurements for inclusive $J/\psi$ production at backward rapidity do not include nuclear effects beyond gluon anti-shadowing.

The ATLAS results on elliptic flow in p+p collisions (Figure 15), the CMS results on $J/\psi$ $v_2$ in p+Pb collisions (Figure 18), and the strong suppression of the excited charmonium (Figure 7) and bottomonium (Figure 18) states at backward rapidity seem to together indicate the formation of quark-gluon plasma in small collision systems. There could be additional effects from comoving particles contributing to quarkonia suppression at backward rapidity, although it seems those effects would be more prominent in the $J/\psi$ nuclear modification results from ALICE (Figure 12).

6. Conclusions

Recent results on quarkonia in small system collisions at RHIC and LHC energies have been presented. The data from PHENIX, STAR, ALICE, and LHCb are consistent with gluon shadowing predictions at forward rapidity, but whether or not additional effects such as energy loss, comover interaction, or dissociation from quark-gluon plasma are present cannot be definitively determined. At backward rapidity, the ALICE $J/\psi$ results and PHENIX open heavy flavor results are consistent with modification only due to gluon anti-shadowing. The PHENIX $J/\psi$ results are consistent with an additional contribution from nuclear absorption, while the LHCb prompt $J/\psi$ measurements show slight suppression from either comoving interactions or dissociation from a hot nuclear medium. Non-zero charm quark elliptic flow was measured by ATLAS and CMS in p+p and p+Pb collisions, respectively. Furthermore, when provided, the transport model predictions can accurately describe quarkonia nuclear modification.

This collection of results suggests quark-gluon plasma formation in small collision systems. However, one of the questions posed in Matsui and Satz’s famous paper is whether or not the $J/\psi$ can “escape from the production region before plasma formation” [96]. Formation time may be an important aspect of small system collision measurements. To quote the preliminary ATLAS results on jet quenching [82]: “It has been proposed that soft (low-momentum) quarks and gluons are only formed on a time scale of 1 fm/c, and thus the high-$p_T$ partons may undergo their virtuality evolution and showering unscathed and fragment in vacuum if the quark–gluon plasma is small, i.e., with a radius < 1–2 fm [97]”.

7. Future Directions

At the time of writing, we are unaware of future upgrades designed for the heavy-ion program at the ATLAS experiment. However, we are aware of significant upgrades coming for the ALICE, CMS, and LHCb heavy-ion programs. The ALICE forward calorimeter (FOCAL), to be installed during the next LHC long shutdown (LS3), offers particle identification and can probe small-x physics [98]. The new CMS minimum–ionising–particle precision timing detector (MTD), part of the Phase II upgrade, will provide particle identification using time-of-flight information [99]. At the LHCb, the fixed target system SMOG II is already installed and recording data for Run 3 [100]. The Magnet Station [101] (for soft particle measurements) and the Mighty Tracker [102] (for central collision measurements) are also scheduled for installation during LS3.

Regarding US experiments, the Electron Ion Collider (EIC) is expected to begin running at Brookhaven National Laboratory around 2032 [103], with a primary goal to study gluon saturation regions at small x. The sPHENIX experiment [104] will come online this spring, with an increased collision rate over PHENIX and the ability to measure the three bottomonium states. In the past thirty years, the heavy-ion community has made exciting discoveries regarding quarkonia production and modification in small and large collision systems. But there is still much to learn and ample opportunity for deeper discovery in the coming decade.