Towards Experimental Confirmation of Quarkonia Melting in Quark–Gluon Plasma: A Review of Recent Measurements of Quarkonia Production in Relativistic Heavy-Ion Collisions

: The dissociation, or “melting”, of heavy quarkonia states due to color charge screening is a predicted signature of quark–gluon plasma (QGP) formation, with a quarkonium state predicted to dissociate when the temperature of the medium is higher than the binding energy of the quarkonium state. A conclusive experimental observation of quarkonium melting coupled with a detailed theoretical understanding of the melting mechanism would enable the use of quarkonia states as temperature probes of the QGP, a long-sought goal in the field of relativistic heavy-ion collisions. However, the interpretation of quarkonia suppression measurements in heavy-ion collisions is complicated by numerous other cold nuclear matter effects that also result in the dissociation of bound quarkonia states. A comprehensive understanding of these cold nuclear matter effects is therefore needed in order to correctly interpret quarkonia production measurements in heavy-ion collisions and to observe the melting of quarkonium states experimentally. In this review, recent measurements of quarkonia production in p A and AA collisions and their state-of-the-art theoretical interpretations will be discussed, as well as the future measurements needed to further the knowledge of cold nuclear matter effects and realize a measurement of quarkonia melting in heavy-ion collisions.


Introduction
Relativistic heavy-ion collisions provide the opportunity to study one of the rarest forms of matter ever created: quark-gluon plasma (QGP).Quarks and gluons, which are normally tightly bound within hadronic bound states of protons and neutrons in nuclei, can be liberated in high-energy collisions of heavy nuclei and result in the formation of the dense partonic plasma phase [1][2][3].QGP is interesting for many reasons in addition to being a novel state of QCD matter.It is predicted to have existed shortly after the Big Bang; therefore, studying it can provide insight into the early formation of the universe, and it also behaves as a nearly perfect hydrodynamic fluid [4,5].Understanding the formation, dynamics, and behavior of QGP would result in significantly improving our understanding of the fundamental physics of the strong nuclear force.
A key property to understand the behavior of any medium is its temperature dependence.The idea to use heavy quarkonia, bound states of charm-anticharm (cc, "charmonia") or beauty-antibeauty (bb, "bottomonia") quarks, as temperature probes in QGP came from a landmark paper by T. Matsui and H. Satz [6].In their paper, Matsui and Satz proposed that quarkonia states would "melt", or be prevented from forming, in a hot deconfined nuclear medium because the free QCD color charges would screen the potential between the two heavy quarks in a mechanism analogous to Debye screening.The screening was predicted to occur when the temperature of the medium exceeded that of the binding energy of the quarkonium state [6].Since different quarkonium states have different binding energies, if one could experimentally observe the melting of quarkonia in a hot nuclear medium, one could infer the temperature of the medium itself.The possibility of experimentally measuring the temperature of QGP has made the observation of the "melting", or dissociation of heavy quarkonia states due to color charge screening, a golden observable in the field of relativistic heavy-ion physics.The idea is that definitive experimental observation of quarkonia melting would enable the use of quarkonium states as temperature probes of quark-gluon plasma.
The straightforward picture of quarkonia melting due to color-charge screening in high-energy nucleus-nucleus ("AA") collisions has evolved over the years as experimental and theoretical studies have revealed that quarkonia production in high-energy collisions is much more complex than initially thought.From the theoretical side, quarkonium production in pp collisions is still not yet completely understood theoretically [7][8][9], and recent theoretical work revealing an imaginary component to the quarkonium potential complicates the picture of the color-charge-screened quarkonium dynamics proposed by Matsui and Satz [10][11][12].From the experimental side, signatures previously thought to only occur in the presence of QGP have been observed even in small collision systems [13], raising questions about whether other QGP effects, including quarkonia melting, could be taking place in small systems as well.Additionally, other effects besides the proposed melting mechanism, commonly called cold nuclear matter (CNM) effects, can cause the dissociation of quarkonia states, and these effects must be fully understood and quantified in order to isolate the signatures of quarkonia dissociation due to melting in QGP [7].Therefore, there are a few requirements that must be met before concluding that quarkonia melting has been experimentally observed.First, it must be proved that no other mechanisms (such as CNM effects) are responsible for quarkonium dissociation.Second, it must be shown that the suppression depends on the binding energy of the measured quarkonia.The first criterion remains much more difficult to satisfy than the second.Disentangling dissociation by melting from dissociation due to other effects remains a challenge.
In this review, recent measurements of quarkonia production in pA and AA collisions will be discussed to highlight the current state of our experimental knowledge of quarkonium suppression in nuclear systems.The photoproduction of quarkonia and the production of quarkonia from beauty hadron decays will not be discussed, as these production mechanisms are different than that for prompt quarkonia and are also not subject to the predicted melting effects.For details on quarkonium production in a vacuum (pp collisions) and comprehensive theoretical introductions, the reader is referred to excellent reviews on these topics [7][8][9][14][15][16].In Section 2, experimental methods for the reconstruction of quarkonia states and the common observables measured to probe quarkonia suppression are briefly reviewed.The discussion of recent experimental results begins in Section 3 with recent results in nucleus-nucleus collisions, where the clearest observations of quarkonium suppression are observed.Then, in Section 4, recent results from proton-nucleus collisions are reviewed, which tend to complicate the interpretation of the results of nucleus-nucleus collisions.

Reconstruction of Quarkonia States
In order to measure the predicted sequential suppression mechanism discussed in the previous section, one must measure both the ground and excited states of the charmonia and bottomonia resonances.For charmonia, the states most experimentally accessible in heavy ion collisions are J/ψ(1S), χ c (1P), and ψ(2S).For bottomonia, the corresponding most accessible states are Υ(1S), Υ(2S), and Υ(3S).With the exception of χ c , the states are typically measured in the dielectron or dimuon decay channel, depending on the experiment.The dilepton decay channels produce clean experimental signals as they produce two charged-particle tracks in a detector that can be precisely measured.Electrons and muons are identified experimentally with dedicated particle identification or muon detectors, respectively.
Fewer experimental measurements exist for the χ c and Υ(3S) states in heavy-ion collisions as these states are particularly challenging to measure experimentally.The χ c state is typically reconstructed in the radiative decay χ cn → J/ψγ, which suffers from significant combinatorial background due to π 0 → γγ decays.The Υ(3S) state is extremely heavy and therefore rarely produced.It is particularly difficult to measure in heavy-ion collisions as the 3S state is predicted to be most affected by dissociation due to cold and hot nuclear matter effects.Therefore, large experimental luminosities are needed to precisely measure Υ(3S) in nuclear collisions.
The production of quarkonia states can be "prompt", where the quarkonium state is produced directly in a hard QCD interaction, or "non-prompt", where the quarkonium is produced in the decay of a beauty hadron.The production of a lower-energy quarkonium state from the decay of a higher-energy quarkonium state is included within the definition of prompt quarkonia; therefore, the only production mechanism for bottomonia states is prompt production.For the charmonia states, however, measuring prompt or nonprompt production probes different physics [7].Prompt charmonia are created directly in a hard QCD interaction in the early stages of a collision, and in the case of heavy-ion collisions are therefore thought to experience the full evolution of the QGP.Non-prompt charmonia are produced in the decays of beauty hadrons, and in this production mechanism, the charmonia are believed to be produced outside the QGP medium [7].In order to probe the dissociation of quarkonia states due to the presence of a nuclear medium, prompt charmonia are the preferred probes.Experimentally, quarkonia can be measured in inclusive, prompt, or non-prompt production.In inclusive production measurements, all candidates for a particular charmonium state in a particular decay channel are measured, and no distinction is made between candidates produced promptly and non-promptly.Prompt and non-prompt charmonia production measurements are performed experimentally by fitting the pseudo-decay time distribution of the charmonia candidates in order to separate candidates produced from b-hadrons, which have larger decay times, from promptly produced charmonia.

Observables for Quarkonia Suppression
Measurements of quarkonia suppression begin by first measuring the cross-section of a quarkonium state in pA, AA, or even pp collisions.Then, one must choose a baseline for which the cross-section measurement will be compared against.Typically, this baseline is either the cross-section of the quarkonium state in pp collisions, or for the case of a higher-resonance quarkonium state, the cross-section of the ground or lower-resonance state measured in the same collision system.
The classic observable for quarkonium suppression is the nuclear modification factor, defined as the ratio of the quarkonium production yield in a nuclear collision system (pA, AA) to the cross-section in pp collisions scaled by the average nuclear overlap, < T AA >, determined from a Glauber model analysis [17]: The nuclear modification factor is typically measured as a function of a variable x, where x is often a kinematic variable, such as the quarkonium transverse momentum p T or rapidity y, or a variable describing the collision centrality, such as the number of participating nucleons in a collision, N part .
Another observable used to probe suppression is the double ratio of quarkonia production in AA collisions relative to that in pp collisions.This ratio is sometimes denoted by ρ: The double ratio of quarkonia production is not as powerful as R AA as it gives only information on the relative suppression between the excited state and the ground state rather than the absolute suppression.However, it is easier to measure experimentally as many factors cancel in the ratio, and combined with an R AA measurement of the ground state, it can be used to derive the excited state R AA .

Bottomonia
The suppression of bottomonia states in PbPb collisions has been studied by ALICE, ATLAS, and CMS [18][19][20][21][22][23][24][25][26].Measurements of elliptic flow of Υ mesons are beyond the scope of this review but are included in the references for the interested reader [27,28].The most precise measurement of bottomonia state suppression in AA collisions comes from CMS, which recently measured all three Υ states in PbPb collisions at √ s = 5.02 TeV [25].The measurements, shown in Figure 1, were performed as a function of centrality and p T and compared to several different theoretical models.A clear sequential suppression is observed for the Υ states, with the Υ(3S) state being the most suppressed, just barely more suppressed than the Υ(2S) state.Both the 2S and 3S states are significantly suppressed compared to the Υ(1S) state.The sequential suppression is observed to be stronger, with the higher Υ states being more suppressed, in central collisions.The suppression measured as a function of p T is constant within the experimental uncertainties [25].
, pp 300 pb The CMS data were compared to several theoretical models, including a coupled Boltzmann model [29], an open quantum system (OQS) model [30], the Heidelberg model [31], the co-mover interaction model (CIM) [32], and the TAMU transport model [33].These models represent a wide variety of approaches to describing the dissociation of quarkonia states in nucleus-nucleus collisions, and their key features will be briefly summarized here.The coupled Boltzmann and TAMU models are both transport models that describe quarkonia production and dissociation in QGP using kinetic rate equations.Typically starting from a space-and momentum-integrated Boltzmann equation, the rate equation contains two transport coefficients that describe the quarkonia dynamics in a nuclear medium [34].The first coefficient is the quarkonia dissociation rate, and the second is the quarkonia yield in the equilibrium limit, which depends on the medium's evolution as a function of temperature.Different transport model calculations vary in their approach to obtaining the two transport coefficients and the medium temperature response.In the TAMU model that is compared to the CMS data in Figure 2, a microscopic model [35] is used to compute the modified quarkonium binding energies in the QGP, and a lattice QCD equation of state is used to compute the medium evolution as a function of temperature [33].Both transport models shown in Figure 2 include contributions from uncorrelated bottomonium recombinations, where a band b-quark pair produced in different initial partonic hard scattering events could recombine to form an Υ meson.The coupled Boltzmann model instead uses two coupled Boltzmann equations: one for an unbound QQ pair (with Q denoting a heavy quark) and one for bound QQ quarkonium states [29].The coupled approach notably provides access to studying the effect of correlated quarkonium regeneration, in which a band b-quark pair produced during the dissociation of a bottomonium state could once again recombine to form a new bottomonium state [29].The open quantum system (OQS) approach treats the quarkonium as a quantum system that interacts with its surrounding environment, in this case the QGP [30,[36][37][38].The reduced density matrix of the quarkonium system is computed, and its equation of motion, the so-called "master equation", is determined [9].The OQS formalism provides a fully quantum description of the quarkonium system, including its dynamics when out of thermal equilibrium with the QGP [30,38] and accounting for quantum correlations between the quarkonium system and the QGP medium [9].The OQS model compared to the CMS data in Figure 2 demonstrates a recent example of combining the OQS framework with an effective field theory that describes the non-relativistic nature of quarkonium, in this case potential non-relativistic QCD (NRQCD) [30,38].The Heidelberg model computes the bottomonium wave functions by solving a temperature-dependent radial stationary Schrodinger equation with an effective potential that includes both real and imaginary components to account for the color screening and damping effects, respectively [31].Finally, the co-mover interaction model describes the dissociation of quarkonia states via their interactions with other "co-moving" particles with similar rapidities in a cold or hot nuclear medium [32].A key parameter in this model is the co-mover density, which is obtained from measurements of the particle multiplicity as a function of rapidity [32].
The double ratio of Υ(3S) to Υ(2S) production yield in PbPb collisions to that in pp collisions provides an excellent testing ground for the theoretical models [25,39].An example of this can be seen in Figure 2, where the measured double ratio and its theoretical predictions are shown as a function of the collision centrality and the quarkonium p T .All models except for the coupled Boltzmann model reasonably describe the data.The OQS predictions are closest to the measured data points, being directly on top of them in several cases.The Heidelberg model also does a good job but overpredicts the ratio at high p T .The TAMU and CIM models are consistently slightly below the data points, and the coupled Boltzmann model is significantly below the data.It would be interesting to see in the future if perhaps adjustments of the relative yields of correlated and uncorrelated bottomonium regeneration in the coupled Boltzmann model could result in better agreement with the CMS data.The R AA of Υ states has also recently been measured at RHIC by STAR [40]. Figure 3 shows the STAR results as a function of centrality compared to results from CMS and to two theoretical models already introduced above, the TAMU transport model [33] and the open quantum system model using potential non-relativistic QCD [36][37][38].It is interesting to observe that Υ(1S) suppression at RHIC is consistent with that observed at the LHC, despite the significantly different center of mass energies.In contrast with Υ(1S), Υ(2S) appears to be slightly more suppressed at RHIC than at the LHC, but the STAR measurement is still consistent with the CMS data within the experimental uncertainties.The theoretical models agree well with the CMS data and the STAR Υ(2S) data; however, both models predict a smaller amount of Υ(1S) suppression at RHIC than what is observed in the data.These recent measurements in the bottomonia sector confirm a sequential suppression pattern but do not yet single out a specific theoretical model that can describe the relative excited state suppression and the suppression dynamics as a function of collision energy.[40] and CMS [22] compared to theoretical models [33,[36][37][38] .See text for details of the theoretical models.Figure by the STAR collaboration, obtained from [40] under American Physical Society copyright and re-used with permission.Image is enlarged with respect to the original version.

Ground State Suppression
The ground state of the charmonia resonances, J/ψ, has been studied extensively in relativistic heavy-ion collisions.Measurements of J/ψ hadroproduction have been performed in PbPb collisions at √ s = 2.76 TeV ( [41][42][43][44][45][46][47]) and 5.02 TeV ( [48][49][50][51][52][53][54][55], in XeXe collisions at √ s = 5.44 TeV ( [56]), in AuAu collisions at √ s = 200 GeV ( [57][58][59][60][61]), and in CuCu collisions at √ s = 200 GeV [58], to name the most recent measurements.Recent flow [62][63][64][65][66][67] and polarization [68] measurements are beyond the scope of this review but are included in the references for the interested reader.Several measurements in non-identical, or AB, nucleus-nucleus collisions have also been recently performed in 3 HeAu collisions by PHENIX [69] and in PbNe collisions by LHCb [70].As previously mentioned in the preceding sections, the discussion here will focus on measurements of prompt charmonia production when available for the aforementioned datasets as prompt measurements are most sensitive to probing dissociation effects due to the nuclear medium.The most recent measurements from RHIC and the LHC, which summarize the current state of the art for J/ψ suppression measurements, will be reviewed and discussed along with some of the most common theoretical approaches used to interpret them.
The most recent measurement of prompt J/ψ production in high-energy nucleusnucleus collisions was performed by ALICE in PbPb collisions at √ s NN = 5.02 TeV [55], where the ALICE results were also compared to measurements by CMS and ATLAS in PbPb collisions at the same collision energy [50,51].The R AA results as a function of J/ψ p T at midrapidity for three different centrality bins are shown in Figure 4.The ATLAS and CMS results span the high-p T region, while the ALICE results cover the low-p T region down to 1.5 GeV/c.All three measurements are consistent in the p T range where they overlap, resulting in a continuous description of the J/ψ R AA from low to high p T .At low p T , the R AA rises dramatically, particularly in the most central collisions, where it is even larger than one.In the next centrality bin from 10-30%, the central values in the two lowest p T bins are slightly lower than in the 0-10% centrality bin but are still consistent with the first centrality bin within the uncertainties.The most noticeable difference is in the p T range between 5 and 10 GeV/c, where the suppression in the 10-30% centrality is slightly less than that seen in the 0-10% bin, consistent with previous measurements that found J/ψ meson suppression to be strongest in the most central collisions.In the final centrality bin, interestingly, J/ψ suppression at low p T is consistent within the uncertainties with the suppression observed by ATLAS and CMS at high p T .
Moving to other collision systems besides PbPb, the inclusive J/ψ R AA has also been recently measured in XeXe collisions by ALICE [56] and in AuAu collisions by STAR [61]. Figure 5 shows the J/ψ R AA measured as a function of centrality in XeXe and PbPb collisions.The XeXe dataset has much larger uncertainties and fewer measurement bins than the PbPb data due to the smaller integrated luminosity available.However, it is still interesting to see that J/ψ suppression in XeXe is largely consistent with that observed in PbPb, despite the fact that the XeXe system is a much lighter collision system (A Xe = 129 and A Pb = 208).The data were compared to the TAMU transport model [71,72], which was previously introduced for the Υ states in Section 3.1.The transport model predictions describe the XeXe data more accurately than the PbPb data, for which they slightly underpredict in the mid-centrality range.In AuAu collisions, the most recent measurement of J/ψ suppression was performed by the STAR collaboration, which measured the J/ψ R AA as a function of p T and centrality [61] .The results shown as a function of p T and binned in five different centrality ranges are shown in Figure 6 and compared to previous results from RHIC and LHC experiments.The J/ψ R AA measured at RHIC is relatively constant as a function of p T within the experimental uncertainties.The difference between the RHIC and LHC results is most prominent at low p T , where the ALICE R AA value is much larger than that measured by STAR.When comparing to the CMS measurements, the RHIC R AA measurement is slightly above the CMS measurement for p T between 8 and 10 GeV/c but is also consistent with the CMS result within the uncertainties.Comparisons with several theoretical models are presented.The Tsinghua [73,74] and TAMU [71,72,75] models are both transport models employing a kinetic rate equation approach to describe charmonium dissociation and regeneration.The transport models tend to agree with each other at low p T , then diverge at high p T , but both are consistent with the data.The remaining collisional dissociation models are only shown in the most central AuAu bin (0-20%), and are based on NRQCD calculations computed with two different values of the charmonium formation time ("t max " and "t min ") [76].The precision of the experimental data cannot yet distinguish between the two formation time scenarios, but this might be possible in the future [77].All of the models also include cold nuclear matter effects on J/ψ production.The enhancement of J/ψ R AA at low p T has been attributed to the production of J/ψ mesons from the recombination of c and c pairs in the extremely hot nuclear medium formed at the LHC [7,9].At high p T , the exact origin of the J/ψ suppression is not yet understood as multiple dissociation mechanisms could be involved, as will be discussed.The picture becomes murkier when moving away from very central collisions towards more peripheral collisions, where the extent of the produced QGP medium and possible interplay with other CNM effects is relatively unknown.The ATLAS and ALICE measurements of the J/ψ R AA were compared to several theoretical predictions, as shown in Figures 7 and 8.The ALICE measurements were compared to three different models: the Statistical Hadronization Model with charm quarks (SHMc) [79], the Boltzmann transport model (BT) [78,80], and a charmonium dissociation model [81,82].The Boltzmann transport model has already been described in the previous section.In the SHMc model, the charm quarks are assumed to thermalize in the QGP medium, and the yields of the charm hadrons are computed with statistical mechanics using a grand canonical partition function [79,83].A key component of the SHMc model is that quarkonia states are only formed during hadronization at the phase transition of the QGP (at the chemical freeze-out temperature T CF ) [83,84].This approach differs from the other models shown in Figure 7, which assume that at least some quarkonia bound states remain intact and can interact with the QGP medium despite the existence of color-screening effects [78,[80][81][82].The dissociation model shown in Figure 7 specifically focuses on describing the dynamics of high p T J/ψs in a QGP and includes J/ψ dissociation contributions from the expected color-charge screening and from collisions of J/ψs with the medium constituents [81].Despite the significantly different physical pictures described by the two models, both the Boltzmann transport model and the SHMc model describe the low p T ALICE data particularly well, while at high values of p T , the dissociation model agrees extremely well with the ALICE data in very central collisions.In semi-central collisions, neither the Boltzmann transport or SHMc models can fully describe the trend of the data.
The ATLAS measurements were compared to color-screening models [81,85] and models involving J/ψ suppression via parton energy loss [86,87].The first color-screening model is the same model that is referred to as the "Dissociation" model in the ALICE plot in Figure 7, which accounts for J/ψ dissociation due to color screening and collisions in the medium [81].The second color-screening model that is compared to the ATLAS data focuses particularly on dissociation due to color exchanges between the QQ state and the QGP medium, which the authors describe as a nuclear absorption-like effect, in addition to the expected melting from color-charge screening [85].Therefore, it is important to clarify that the "color-screening" models compared to the ATLAS data do not only include quarkonia dissocation due to color screening, but they also include dissociation due to non-melting effects as well.The first energy loss model referred to in the ATLAS figure describes quarkonia suppression in AA collisions as being predominantly due to the radiative energy loss of the charmonia in the QGP [86].The radiative energy loss model is based on the idea that quarkonia suppression in QGP could have the same underlying origin as the phenomenon of jet "quenching", in which jets, high-energy sprays of particles produced from the hadronization of a single quark or gluon, are suppressed or "quenched" in the QGP compared to pp collisions [86].Phenomenologically, the quenching (or suppression) depends on the original p T of the quark or gluon (or charmonium state) and an effective color factor that accounts for the different radiation patterns between quarks and gluons [86].It is indeed interesting that this quenching model gives almost identical predictions as the first color-screening model despite the very different underlying theoretical mechanisms, and both are in very good agreement with the experimental data.The second energy loss model compared to the ATLAS data shares the idea that radiative energy loss is the dominant mechanism explaining J/ψ suppression, but rather than comparing J/ψ suppression with the quenching of jets, in this approach, it is instead compared with the suppression of individual hadrons [87].It is also found in this model that J/ψ suppression and hadron suppression can be described by the same underlying mechanism.One caveat to the energy loss models is that they are expected to be most relevant at high p T (>≈ 10 GeV/c); therefore, they are not expected to describe the observed J/ψ suppression at low p T where melting could be playing an important role.However, it is certainly interestingly to observe that the ATLAS data at high p T are consistent with both dissociation due to color screening and suppression due to energy loss.No single model in Figures 7 and 8
The most recent measurement comparing the suppression of ψ(2S) to J/ψ mesons was performed by the ALICE collaboration at √ s NN = 5.02 TeV [88].The suppression factor R AA was measured as a function of p T and centrality and compared to CMS results measured at the same center of mass energy [51].Figure 9 shows the R AA of J/ψ and ψ(2S) mesons measured by ALICE and CMS as a function of the quarkonia transverse momentum p T .Note that the ALICE results are for inclusive J/ψ and ψ(2S), while the CMS results are for prompt J/ψ and ψ(2S).The ALICE and CMS results nicely complement each other, with the ALICE results covering the extremely low p T region and the CMS results spanning the high p T region up to p T = 30 GeV/c.The data clearly show that ψ(2S) has a lower R AA across the entire p T range measured compared to J/ψ, indicating that ψ(2S) is indeed more suppressed than J/ψ in PbPb collisions.The data were compared to theoretical predictions based on the TAMU transport model [72], which was previously introduced in Section 3.1.
The model calculations describe the data quite well, in particular accurately capturing the rise at low p T attributed to charmonium recombination.The R AA of prompt J/ψ and ψ(2S) has also been measured as a function of the collision centrality by both ALICE and CMS [51,88].The measured R AA distributions are shown in Figure 10, with the CMS measurement shown on the left and the ALICE measurement shown on the right.When comparing the two measurements as a function of centrality, it is important to note that the experiments probe different charmonia p T regions.This explains why the two measurements are consistent with each other when integrated over centrality and plotted as a function of p T but appear to disagree when integrated over p T and plotted as a function of centrality.Nevertheless, although the absolute value of R AA for a given value of N part cannot be directly compared between the two measurements, the trends are still interesting.The ALICE measurement, which has a lower charmonia p T range than the CMS measurement, flattens out towards the most central (highest < N part >) collisions.The CMS results, however, show a generally decreasing trend of R AA towards the most central collisions for prompt J/ψs.The ψ(2S) results in the CMS measurement have larger uncertainties, so the trend with centrality is not as clear, but they could hint at a possible decrease with centrality as well.[88], the ALICE J/ψ data is from [52], and the CMS data are from [51].Comparisons are shown to predictions using the TAMU transport model [72].[88], with comparisons to theoretical models [72,79,89].The ALICE J/ψ data is from [49].The double ratio of ψ(2S) to J/ψ yields in PbPb collisions to the ratio in pp collisions has been measured by ALICE, ATLAS, and CMS at √ s = 5.02 TeV [48,50,88].The ALICE measurements of the double ratio as a function of p T indicate that ψ(2S) is suppressed by roughly a factor of 40-50% in PbPb collisions compared to pp collisions [88].The CMS measurement hints at a slightly larger suppression at midrapidity, but in the most forward rapidity bin probed, 1.6 < |y| < 2.4, the suppression is consistent with that observed by ALICE [48,88].ATLAS measured the double ratio as a function of centrality for both prompt and non-prompt charmonia.The double ratio for the non-prompt charmonia was observed to be consistent with one, supporting the interpretation that the non-prompt charmonia, which originate from b-hadron decays, hadronize outside the medium and are therefore unaffected by the suppression mechanisms responsible for dissociating prompt charmonia in the medium [50].

ALI-PUB-528412
A key missing component in the picture of sequential charmonium suppression in AA collisions is the χ c (1P) state.With a binding energy in between that of J/ψ and ψ(2S), a measurement of its suppression is needed to fully confirm the sequential suppression of charmonia states.As with the Υ states, no single model can fully describe all of the experimental data on the charmonia, which likely indicates, as the variety of available models also suggests, that multiple mechanisms contribute simultaneously to charmonia dissociation in QGP.

Recent Results from pA Collisions 4.1. Bottomonia
Bottomonia production in pA collisions has been predominantly studied at the LHC, with measurements reported by all LHC experiments in pPb collisions [90][91][92][93][94][95].The most recent measurement of the nuclear modification of bottomonia states at the LHC was performed by CMS [95]. Figure 11 shows the R pPb for the Υ(1S), Υ(2S), and Υ(3S) states measured as a function of rapidity for p T < 30 GeV/c.A sequential suppression of the upsilon states is observed, with Υ(3S) being more suppressed than Υ(2S), which is more suppressed than Υ(1S).The sequential suppression is flat within the experimental uncertainties as a function of the quarkonia p T and y.The data in Figure 11 were compared to predictions from the co-mover interaction model (CIM) [32] using two different sets of nuclear parton distribution functions (nPDFs), nCTEQ15 [96], and EPS09 at leading order (LO) [97].The CIM was previously introduced in Section 3.1 and describes quarkonia dissociation via collisions with co-moving particles at similar rapidities as the quarkonium [32].The use of different nPDFs is to test different descriptions of the initial nuclear state.Obtained from fits to a wide variety of data, nPDFs quantify the difference in the partonic structure of a bound nucleus compared to a hypothesis of the partonic structure in A non-interacting protons and neutrons, where A is the atomic mass number.If a bound nucleus were simply a superposition of non-interacting protons and neutrons, the nPDFs would not show significant modification relative to the A-scaled proton PDFs.However, indeed, studies of vast amounts of experimental data have revealed that nuclei are not simple superpositions of non-interacting hadrons; they have their own complex structure, of which little is understood on a partonic level [96].One such example of this structure is the existence of nuclear "shadowing", in which the nuclear PDF is observed to be suppressed in certain partonic kinematic regions compared to the same regions in a proton PDF for reasons not yet fully understood.Recent progress on nPDF interpretation and fitting is described in the excellent review in [98].Both nPDF sets used in conjunction with the CIM model are fitted to deep inelastic scattering and Drell-Yan data measured with nuclear targets and also include inclusive π 0 production data measured in dAu collisions at RHIC [96,97].The nCTEQ15 nPDFs are fitted to fewer data points than the EPS09 nPDFs, and the determination of the uncertainty approaches also varies between the two collaborations [96].The CIM model with either nPDF set does a good job of describing the CMS Υ data, with perhaps the only exception being a slight discrepancy with the Υ(3S) data at negative rapidity.
The Υ(1S) R pPb data were also compared to predictions involving an energy loss mechanism [99] as well as predictions that only included nuclear shadowing effects using the EPS09 nPDFs computed at next-to-leading-order [97,100] shown in Figure 12.The energy loss model depicted here specifically accounts for the fully coherent energy loss (FCEL) mechanism that occurs due to the interference of gluon emission amplitudes before and after a partonic hard scattering event [101,102].The data are consistent with all of the compared theory models, indicating that the observed Υ suppression is consistent with being produced as a result of CNM effects, but the exact effects responsible for the suppression remain unknown.[95] as a function of y compared to predictions from the co-mover interaction model (CIM) [32] with nCTEQ15 [96] and EPS09 LO [97] nPDFs.Figures are by the CMS Collaboration, obtained from [95] under CERN copyright, and licensed under CC BY 4.0.The size of the figures is reduced with respect to the original version.[95] compared to (left) predictions involving nuclear shadowing effects [100] and (right) predictions involving shadowing and energy loss effects [99].Both predictions use the EPS09 NLO nPDFs [97].Figures are by the CMS Collaboration, obtained from [95] under CERN copyright, and licensed under CC BY 4.0.The size of the figures is reduced with respect to the original version.
Recent measurements of Υ production in √ s = 8.16 TeV pPb collisions were also performed by the ALICE experiment [94].Figure 13 shows the ALICE and LHCb measurements of the Υ(1S) R pPb compared to several theoretical models.The first is a model using the EPS09NLO nPDFs to describe the nuclear shadowing effects combined with the color evaporation model (CEM) to describe Υ production [100].In the CEM, the mass of a QQ pair is restricted to be less than the open heavy flavor meson threshold, and the total cross-section for a quarkonium state given this mass constraint is computed from a fraction F of the total QQ pairs satisfying this criteria [100].The next models were previously introduced during the discussion of Figure 12 and include fully coherent energy loss model predictions [103] with and without the additional inclusion of nuclear shadowing effects with the EPS09 nPDFs [104].The EPPS16 reweighted and nCTEQ15 reweighted predictions are, respectively, the EPPS16 [105] and nCTEQ15 [96] nPDF sets reweighted with heavy flavor hadron production data from the LHC to study the effect of the heavy flavor data on constraining the gluon nuclear PDF [106].The final model shown is the co-mover model with nCTEQ15 nPDFs for shadowing contributions [32].As seen with the CMS results, the data from ALICE and LHCb are consistent with models involving several different CNM effects.The co-mover model is particularly promising in that it can describe the suppression of all three Υ states, and it appears to do a slightly better job of describing the data at negative rapidity than the other models.

Charmonia 4.2.1. Ground State Suppression
Following the strong suppression of J/ψ mesons observed in AA collisions, it is natural to ask what we currently know about J/ψ suppression in pA collisions.At the LHC, measurements of J/ψ R pPb have been performed at √ s = 5.02 TeV and 8.16 TeV [93,[111][112][113][114][115][116][117][118][119][120].As with the AA results, measurements of prompt J/ψ production will be highlighted here as these are most sensitive to probing J/ψ dissociation mechanisms.The LHC measurements will be briefly introduced, and then their theoretical comparisons will be discussed since the data are compared to many of the same theoretical models.Figure 14 shows the most recent prompt J/ψ R pPb results from ALICE [120] compared to previous LHC measurements at midrapidity, as well as several theory predictions and plotted as a function of J/ψ p T .While one rapidity bin of the CMS measurement is shown in Figure 14, the measurement was also performed in multiple rapidity bins to probe the nuclear effects on J/ψ production as a function of both p T and rapidity, as shown in Figure 15.Due to its forward kinematic acceptance, LHCb has measured the prompt J/ψ R pPb in a complementary rapidity range to that of the midrapidity LHC experiments, probing more positive and negative rapidity values in the center-of-mass frame [112].The LHCb R pPb results are shown as a function of p T for two different rapidity bins in Figure 16.The LHCb and ALICE R pPb measurements are shown as a function of the center of mass rapidity in Figure 17.

ALI-PUB-561246
Figure 14.The prompt J/ψ R pPb measured by ALICE [118,120], ATLAS [121], and CMS [111] as a function of the J/ψ p T and compared to several theoretical predictions [104][105][106][107]122,123].See text for descriptions of the theory models.All of the R pPb results in Figures 14-17 are compared to theoretical predictions, most of which incorporate different cold nuclear matter effects.Several of the nuclear PDFs used in the models in Figures 14-17 have already been previously introduced, such as the EPS09 at LO and NLO [105] and the nCTEQ15 [96] and EPPS16 [105] nPDFs with reweighting to include the impact of LHC heavy flavor data [106,107].A new addition is the nCTEQ15HQ set [122], which included open heavy flavor and quarkonium measurements in pPb data from the LHC in a global nPDF analysis (the "HQ" in the name stands for heavy quark).The HELAC-Onia models compared to the LHCb data in Figure 16 are based on the HELAC-Onia package, which computes the matrix elements involved in collinear NRQCD factorisation predictions of quarkonia production [108,109].A color glass condensate (CGC) approach is also introduced as a model comparison to the J/ψ R pPb data [124][125][126].The CGC model is an effective field theory applicable when the saturation scale, the scale at which the gluon density in a nucleus "saturates" and reaches its maximum limit, is much larger than the energy scale associated with nonperturbative QCD, denoted by Λ QCD [124].The cc pair is produced from the radiation of a "probe" parton, and the pair proceeds to interact with the dense gluonic system in the target nucleus [125,127].In contrast to the rest of the models, the transport model predictions [123] in Figure 14 involve a hot nuclear medium and the treatment of quarkonia production and dissociation using a rate equation approach as described in Section 3.1.

ALI-PUB-530723
Figure 17.The J/ψ R pPb measured by ALICE [118] and LHCb [119] as a function of rapidity and compared to several theoretical models [96,97,[105][106][107][108][109]124,125,128,129].See text for details of the models.Several interesting features can be noted from comparing the midrapidity results in Figures 14 and 15 with the forward and backward rapidity results in Figure 16.At midrapidity and high p T , the R pPb is slightly larger than one.This can be most clearly seen in the CMS results shown in Figure 15, in which the R pPb is observed to be larger than one in most of the measured p T bins in −0.9 < y CM < −1.93.It is when approaching the more forward rapidities that the value of R pPb begins to drop below one, consistent with what is observed in the LHCb data.The R pPb is smaller than one for almost the entire p T range probed in the forward region and for p T < 3 GeV at midrapidity.The R pPb plotted as a function of rapidity shows a decreasing trend from negative to positive rapidity, as shown in Figure 17.Another observation is that the data cannot be described by a single nuclear PDF set, as can be seen most clearly with the LHCb data.The data are consistent with the models, which seem to suggest that nPDF modifications, at a minimum, are needed to describe the data but that other CNM effects, such as energy loss or saturation effects, could also be at play.
RHIC data are also extremely valuable for studying cold nuclear matter effects as quarkonia measurements in many nuclear systems, collision energies, and regions of phase space are needed to begin to disentangle the different effects and understand their underlying mechanisms.An example of this can be seen in the recent measurement from PHENIX of the J/ψ R pAl , R pAu , and R HeAu [69].The R AB results from PHENIX are shown for each collision system as a function of p T in Figure 18 and as a function of rapidity in Figure 19.The theoretical models, which have already been discussed previously, include different parameterizations of nuclear shadowing effects via the EPPS16 nPDFs [100,105], the EPS09 nPDFs [105], and the reweighted EPPS16 and nCTEQ15 nPDFs [106,107] and a transport model that includes hot medium effects on charmonia production [75,123].Contrary to the LHC results, the PHENIX results shown here begin to differentiate between theoretical models, as can be seen for the observed values of the R AB at negative rapidity in p+Au and 3 He+Au, which do not agree with the predictions from Vogt and Shao.These predictions did not include nuclear absorption effects, which is expected to be significant at RHIC energies [130].When adding a nuclear absorption model [131] to the same calculations (predictions labeled with "+Abs" in the figure), the agreement between the data and predictions is greatly improved.This study presents a nice example of how a specific cold nuclear matter effect can begin to be isolated, or at least identified as the dominant effect, in different nuclear systems in order to better understand its underlying mechanism.The STAR collaboration has also recently measured the J/ψ R pAu , as can be seen in Figure 20 [132].The STAR pAu results are compared to the J/ψ R AA measurements by STAR in AuAu collisions, to dAu results measured by PHENIX, and to several theoretical calculations.The pAu and dAu results are consistent with each other over the entire measured p T range, which suggests that the measured nuclear modification is primarily due to nuclear effects associated with the Au nucleus.There are several similarities between the R pAu measured at midrapidity and the R pPb measured at midrapidity.For both cases, the nuclear modification factor is consistent with one starting around 3 GeV/c.The small enhancement above one at high p T that was observed in the ATLAS and CMS data is not observed in the ALICE or STAR data, possibly due to the larger statistical uncertainty for the latter results.The theoretical comparisons include effects due to nuclear shadowing [97,[106][107][108][109], a CGC model [133,134], a transport model including a hot nuclear medium effect [72,123], an energy loss model [128], and a co-mover model [135].Notably, the STAR data are consistent with predictions involving cold nuclear matter effects and the TAMU transport model, which involves the formation of a QGP medium.The co-mover model predictions seem to do particularly well at describing the data at low p T , but the other models do a better job of describing the data over the entire measured p T range.It is worth mentioning that particularly for the J/ψ results discussed in this section, the data can be described by many models that describe individual CNM effects, usually in combination with nuclear shadowing.However, it is likely the case that multiple CNM effects occurring simultaneously are responsible for the observed J/ψ dissociation.Understanding the relative contributions of various CNM effects to the observed suppression remains a challenging but important goal.The nuclear modification factor R pAu for J/ψ mesons measured as a function of p T by STAR [61,132] and PHENIX [136] and compared to several theoretical models [72,96,[106][107][108][109]123,128,[133][134][135].

Excited State Suppression
Proton-nucleus systems are ideal for studying excited-state charmonium production as these states are generally statistically limited in nucleus-nucleus collisions.Recent measurements in pA collisions have studied χ c production, which is particularly interesting as it has never been studied in AA collisions.As previously discussed in Section 2, the prompt J/ψ yield measured experimentally includes J/ψ production via the decays of higher charmonium states, including χ c .LHCb measured for the first time in pPb collisions at the LHC the fraction of J/ψ mesons produced from χ c decays [137], where the total χ c yield included the χ c1 and χ c2 states.The result is shown in Figure 21 as a function of J/ψ p T .The pPb measurement is consistent with the fraction measured in pp collisions except in the lowest p T bins at negative rapidity; however, the increased fraction in these bins is compatible with the observed ψ(2S) suppression in pPb collisions [137].The production of ψ(2S) mesons in pPb collisions has been measured by all four LHC experiments [93,117,[139][140][141][142][143][144]. Figure 22 shows the ψ(2S) R pPb measured by CMS as a function of y in three different p T bins [142].The ψ(2S) suppression is most evident in the lowest p T bin probed from 4 to 6.5 GeV/c, where the ψ(2S) is significantly suppressed with respect to the J/ψ, particularly at negative rapidity.In the middle p T bin, from 6.5 to 10 GeV/c, the ψ(2S) appears significantly suppressed at negative rapidity but is consistent with J/ψ and no suppression for y > 1.In the highest p T bin probed, the ψ(2S) suppression is consistent with one, although it is systematically lower than the J/ψ R pPb values.It is interesting that the R pPb values for J/ψ are consistently slightly larger than one, while the values for ψ(2S) tend to be smaller than one except in the highest p T bin.
LHCb has also measured the prompt ψ(2S) R pPb and compared it to several theoretical models [140].The LHCb measurement as a function of rapidity is shown in Figure 23.The LHCb results clearly show that the ψ(2S) is more suppressed than the J/ψ at negative rapidity than at positive rapidity, which is also consistent with the trend observed by CMS in the p T bins where the measurements overlap.The LHCb data are compared to models involving shadowing via the EPS09 PDFs at LO [129,145] and NLO [146] and to energy loss predictions with and without nuclear shadowing effects included [103].Interestingly, while J/ψ suppression can be described by the theoretical models, none of them successfully describe ψ(2S) suppression.The same difficulty of theoretical models in describing both J/ψ and ψ(2S) suppression was observed by ALICE when they measured the inclusive ψ(2S) R pPb [143].Figure 24a shows the ALICE measurements compared to several theoretical models, employing some of the same cold nuclear matter effects as were included in the models compared to the LHCb data.The models include shadowing effects [106,146], CGC models [126,134], an energy loss model [99], and the co-mover model [135].The co-mover model compared to the data in Figure 24b appears to be the only model able to describe both J/ψ and ψ(2S) suppression, at least calculating the approximate order of magnitude for the ψ(2S) suppression correctly.The nuclear modification factor R pPb for prompt J/ψ [119] and ψ(2S) [140] mesons measured as a function of rapidity and compared to several theoretical models [97,103,129,145,146].The nuclear modification factor of inclusive J/ψ [113] and ψ(2S) [143] mesons as a function of rapidity compared to several theoretical models.(a) depicts comparisons to models using different nPDFs [96,97,104,106] as well as energy loss [99,104] and CGC [126] models and (b) depicts comparisons to CGC [134] and co-mover [135]  Recent studies of ψ(2S) suppression have also been performed at RHIC by PHENIX in pAl and pAu collisions [147].The PHENIX results also show significant ψ(2S) suppression at negative rapidity that cannot be explained by models using different nuclear PDFs, specifically the reweighted EPPS16 and nCTEQ15 nPDFS [106,107], as shown in Figure 25.The authors note that nuclear absorption cannot explain the observed suppression as the initial-state CNM effects should be similar for J/ψ and ψ(2S).It would, however, still be interesting to confirm this by comparing to models including a nuclear absorption component, as was conducted in [69].It would also be interesting to see if the co-mover model can reproduce the observed ψ(2S) suppression in the backward region as it appears to be able to do for the ALICE data shown in Figure 24, which could also not be described by shadowing alone.The PHENIX nuclear modification factor measurements were also compared to recent measurements from ALICE [117] as a function of centrality and transport model predictions [123], which are shown in Figure 26.Note that in Figure 26, the ALICE nuclear modification factor as a function of centrality is denoted by Q pPb due to potential biases in the determination of centrality in pPb collisions.The transport model predictions describe the data better in the forward rapidity region than in the backward rapidity region, although the backward predictions are still consistent with the data points within the large uncertainties.It is important to note here that the ALICE data shown in Figure 26 were also well described by co-mover model predictions, again showing that the observed suppression is consistent with predictions involving hot or cold nuclear matter effects.

Conclusions
Significant progress has been made both experimentally and theoretically to better understand quarkonium production and modification in nuclear systems.The measurements presented here show that the production of quarkonia in nuclear systems is indeed suppressed with respect to pp collisions, and the suppression in nucleus-nucleus collisions where QGP is expected to be present is indeed stronger than that observed in smaller proton-nucleus systems.There is experimental evidence of the sequential suppression of quarkonia states as a function of their binding energy for both the bottomonia and charmonia systems.However, the experimental measurements and theoretical descriptions cannot yet definitively conclude that the observed suppression in AA collisions is due to quarkonium melting.Other effects, such as parton energy loss and co-mover interactions, can also describe the data within the current experimental uncertainties.In order to better understand the complex quarkonium dynamics in AA collisions, the cold nuclear matter effects that modify quarkonium production in pA collisions, and the quarkonium production mechanisms in pp collisions, the baseline for comparisons to AA collisions, must be more thoroughly understood.Future measurements in a variety of collision systems combined with the increased luminosities expected with future RHIC and LHC data will enable precise and differential measurements that will help to improve our knowledge of quarkonia production and dissociation in nuclear systems.

Figure 1 .
Figure 1.The R AA of the Υ(1S), Υ(2S), and Υ(3S) states measured by CMS as a function of (left) centrality and (right) quarkonia p T .Figures are by the CMS Collaboration, obtained from [25] under CERN copyright, and licensed under CC BY 4.0.Images are cropped from the original versions.

Figure 2 .
Figure 2. The R AA of the Υ(1S), Υ(2S), and Υ(3S) states measured by CMS compared to several theoretical models as a function of (a) centrality and (b) quarkonia p T .See text for details of the theoretical models.Figures by the CMS Collaboration, obtained from [39] under CERN copyright.The size of the images is reduced with respect to the original versions.

Figure 3 .
Figure 3.The R AA of Υ(1S) (top) and Υ(2S) (bottom) mesons measured by STAR[40] and CMS[22] compared to theoretical models[33,[36][37][38] .See text for details of the theoretical models.Figure by the STAR collaboration, obtained from[40] under American Physical Society copyright and re-used with permission.Image is enlarged with respect to the original version.

Figure 5 .
Figure5.The R AA for J/ψ mesons in XeXe collisions[56] compared to PbPb collisions[49], measured as a function of centrality and compared to transport model predictions[71,72].Figure by the ALICE Collaboration, obtained from[56] under CERN copyright, and licensed under CC BY 4.0.The size of the image is reduced with respect to the original version.

Figure 6 .
Figure 6.The J/ψ R AA measured by STAR [61] as a function of p T and binned in centrality.AuAu results from RHIC [57-59] are compared with PbPb results from the LHC [46,63].The results are compared to several theoretical predictions [73-76,78].Figure by the STAR Collaboration, obtained from [61], and licensed under CC BY 4.0.The figure is reduced with respect to the original version.
can describe the suppression over the full measurement range.

Figure 7 .
Figure 7.The J/ψ R AA measured by ALICE [55] as a function of p T and compared to predictions from several theoretical models [78-82].The data and theory comparisons are shown for two centrality bins: (a) 0-10% and (b) 30-50%.Figures are by the ALICE Collaboration, obtained from [55] under CERN copyright, and licensed under CC BY 4.0.The size of the images is reduced with respect to the original versions.

Figure 8 .
Figure 8.The prompt J/ψ R AA measured by ATLAS [50] as a function of p T and compared to predictions from theoretical models [81,85-87].Figure by the ATLAS Collaboration, obtained from [50] under CERN copyright, and licensed under CC BY 4.0.The figure size is reduced with respect to the original version.

Figure 9 .
Figure 9.The nuclear modification factor R AA for J/ψ and ψ(2S) mesons in PbPb collisions at √ s = 5.02 TeV as a function of the quarkonium transverse momentum p T .The ALICE data are inclusive mesons, while the CMS data are prompt mesons.The figure and ALICE ψ(2S) data are from[88], the ALICE J/ψ data is from[52], and the CMS data are from[51].Comparisons are shown to predictions using the TAMU transport model[72].Figure by the ALICE Collaboration, obtained from [88] under CERN copyright, and licensed under CC BY 4.0.Figure is enlarged with respect to original version.

Figure 10 .
Figure 9.The nuclear modification factor R AA for J/ψ and ψ(2S) mesons in PbPb collisions at √ s = 5.02 TeV as a function of the quarkonium transverse momentum p T .The ALICE data are inclusive mesons, while the CMS data are prompt mesons.The figure and ALICE ψ(2S) data are from[88], the ALICE J/ψ data is from[52], and the CMS data are from[51].Comparisons are shown to predictions using the TAMU transport model[72].Figure by the ALICE Collaboration, obtained from [88] under CERN copyright, and licensed under CC BY 4.0.Figure is enlarged with respect to original version.
Figure 10.The nuclear modification factor R AA for prompt J/ψ and ψ(2S) mesons in PbPb collisions at √ s = 5.02 TeV as a function of the centrality (here defined with the average number of particles participating in the interaction, N part ).(a) Measurement by CMS, from [51]. Figure by the CMS Collaboration, obtained from [51] under CERN copyright, and licensed under CC BY 4.0.The figure size is reduced with respect to the original version.(b) Measurement by ALICE, from[88], with comparisons to theoretical models[72,79,89].The ALICE J/ψ data is from[49].Figure by the ALICE Collaboration, obtained from[88] under CERN copyright, and licensed under CC BY 4.0.The figure size is reduced with respect to the original version.Note that the color scheme for J/ψ and ψ(2S) is switched between the two figures.

Figure 11 .
Figure 11.The nuclear modification factors R pPb measured for Υ(1S), Υ(2S), and Υ(3S) mesons[95] as a function of y compared to predictions from the co-mover interaction model (CIM)[32] with nCTEQ15[96] and EPS09 LO[97] nPDFs.Figures are by the CMS Collaboration, obtained from[95] under CERN copyright, and licensed under CC BY 4.0.The size of the figures is reduced with respect to the original version.

Figure 12 .
Figure12.The nuclear modification factor R pPb measured for Υ(1S) mesons[95] compared to (left) predictions involving nuclear shadowing effects[100] and (right) predictions involving shadowing and energy loss effects[99].Both predictions use the EPS09 NLO nPDFs[97].Figures are by the CMS Collaboration, obtained from[95] under CERN copyright, and licensed under CC BY 4.0.The size of the figures is reduced with respect to the original version.

Figure 15 .
Figure 15.The J/ψ R pPb measured by CMS [111] as a function of p T and y.See text for descriptions of the theory models [96,97,100,107-109]. Figure by the CMS Collaboration, obtained from [111] under CERN copyright, and licensed under CC BY 4.0.The size of the figure is reduced with respect to the original version.

Figure 16 .
Figure 16.The prompt J/ψ R pPb measured by LHCb [112] as a function of the rapidity in the center of mass frame y * in two different rapidity bins corresponding to (a) the pPb beam configuration with the proton beam entering the detector and (b) the Pbp beam configuration with the Pb beam entering the detector.See text for details of the theoretical comparisons [108,109,125,126].Figures by the LHCb Collaboration, obtained from [112] under CERN copyright, and licensed under CC BY 4.0.The size of the figures is reduced with respect to the original version.
Figure 17.The J/ψ R pPb measured by ALICE[118] and LHCb[119] as a function of rapidity and compared to several theoretical models[96,97,[105][106][107][108][109]124,125,128,129].See text for details of the models.Figure by the ALICE Collaboration, obtained from [118] under CERN copyright, and licensed under CC BY 4.0.The figure size is enlarged with respect to the original version.

Figure 18 .
Figure 18.The J/ψ nuclear modification factor as a function of p T measured by PHENIX [69] in three different collision systems: pAl (a) , pAu (b) , and 3 HeAu (c) .The top row shows the results at negative rapidity, while the bottom row shows the results at positive rapidity.See text for details of the theoretical predictions [96,97,100,105-109,123] .Figures by the PHENIX Collaboration, obtained from [69], and licensed under CC BY 4.0.The size of the figures is reduced with respect to the original version.

Figure 19 .Figure 20 .
Figure 19.The J/ψ nuclear modification factor as a function of y measured by PHENIX [69] in three different collision systems: pAl (a), pAu (b), and 3 HeAu (c).See text for details of the theoretical predictions [96,100,105-109,130,131]. Figure by the PHENIX Collaboration, obtained from [69], and licensed under CC BY 4.0.The size of the figure is reduced with respect to the original version.

Figure
Figure by the STAR Collaboration, obtained from [132], and licensed under CC BY 4.0.The size of the figure is reduced with respect to the original version.

Figure 21 .
Figure 21.The fraction of prompt J/ψ mesons produced from χ c decays measured in pp [138], pPb and Pbp collisions [137] as a function of J/ψ p T .Figure by the LHCb Collaboration, obtained from [137] under CERN copyright, and licensed under CC BY 4.0.The size of the figure is reduced with respect to the original version.

Figure 22 .Figure 23 .
Figure 22.The nuclear modification factor R pPb for prompt J/ψ [111] and ψ(2S) [142] mesons measured as a function of rapidity and p T .The top left figure shows the results in the lowest p T bin probed, 4-6.5 GeV/c; the top right figure shows the results in p T = 6.5-10GeV/c; and the bottom figure shows the results in p T = 10-30 GeV/c.Figures by the CMS Collaboration, obtained from [142] under CERN copyright, and licensed under CC BY 4.0.The size of the figures is reduced with respect to the original version.
Figure23.The nuclear modification factor R pPb for prompt J/ψ[119] and ψ(2S)[140] mesons measured as a function of rapidity and compared to several theoretical models[97,103,129,145,146]. Figure by the LHCb Collaboration, obtained from [140] under CERN copyright, and licensed under CC BY 4.0.The size of the figure is enlarged with respect to the original version.

Figure 24 .
Figure24.The nuclear modification factor of inclusive J/ψ[113] and ψ(2S)[143] mesons as a function of rapidity compared to several theoretical models.(a) depicts comparisons to models using different nPDFs[96,97,104,106] as well as energy loss[99,104] and CGC[126] models and (b) depicts comparisons to CGC[134] and co-mover[135] model predictions.See text for details of the theoretical models.Figures by the ALICE Collaboration, obtained from [143] under CERN copyright, and licensed under CC BY 4.0.The size of the figures is reduced with respect to the original version.
Figure24.The nuclear modification factor of inclusive J/ψ[113] and ψ(2S)[143] mesons as a function of rapidity compared to several theoretical models.(a) depicts comparisons to models using different nPDFs[96,97,104,106] as well as energy loss[99,104] and CGC[126] models and (b) depicts comparisons to CGC[134] and co-mover[135] model predictions.See text for details of the theoretical models.Figures by the ALICE Collaboration, obtained from [143] under CERN copyright, and licensed under CC BY 4.0.The size of the figures is reduced with respect to the original version.

Figure 25 .
Figure 25.The J/ψ and ψ(2S) nuclear modification factors measured as a function of rapidity by PHENIX in (a) pAl and (b) pAu collisions and compared to theoretical predictions [96,105-109].

Figures
Figures by the PHENIX collaboration, obtained from[147] under American Physical Society copyright, and re-used with permission.The size of the figures is reduced with respect to the original version.

Figure 26 .
Figure 26.The ψ(2S) nuclear modification factor measured as a function of centrality by PHENIX (R pAu )[147] and ALICE (Q pPb )[117] and compared to transport model predictions[123].Figure by the PHENIX collaboration, obtained from[147] under American Physical Society copyright, and re-used with permission.The size of the figure is reduced with respect to the original version.