In-Medium Charmonium Production in Proton-Nucleus Collisions

We study charmonium production in proton-nucleus ($p$-A) collisions focusing on final-state effects caused by the formation of an expanding medium. Toward this end, we utilize a rate equation approach within a fireball model as previously employed for a wide range of heavy-ion collisions, adapted to the small systems in $p$-A collisions. The initial geometry of the fireball is taken from a Monte-Carlo event generator where initial anisotropies are caused by fluctuations. We calculate the centrality and transverse-momentum dependent nuclear modification factor ($R_{p{\rm A}}$) as well as elliptic flow ($v_2$) for both $J/\psi$ and $\psi(2S)$ and compare them to experimental data from RHIC and the LHC. While the $R_{p{\rm A}}$s show an overall fair agreement with most of the data, the large $v_2$ values observed in $p$-Pb collisions at the LHC cannot be accounted for in our approach. While the former finding generally supports the formation of a near thermalized QCD medium in small systems, the discrepancy in the $v_2$ suggests that its large observed values are unlikely to be due to the final-state collectivity of the fireball alone.


Introduction
Quarkonia, bound states of a heavy quark (Q=b, c) and its antiquark, are pristine objects to study the properties of the fundamental color force. The discovery of charmonia and bottomonia in the 1970s led to the development of the Cornell potential [1] which remains valid to date. To study its modifications in QCD media, large-scale experimental and theoretical efforts are ongoing to measure and interpret the production systematics of quarkonia in ultra-relativistic heavy-ion collisions (URHICs) [2][3][4][5]. An interplay of suppression processes in a quark-gluon plasma (QGP) and subsequent regeneration in the later stages of the fireball evolution is, in principle, required to arrive at a comprehensive description of quarkonium observables in URHICs. In particular, the significance of regeneration reactions has been validated by an excitation function for the J/ψ which transits from a regime of strong suppression at SPS and RHIC energies to significantly less suppression at the LHC, with the extra yield mostly concentrated at low transverse momenta (p T ) (cf. Ref. [6] for a recent overview). However, when accounting for the so-called coldnuclear matter (CNM) effects, which, in particular, cause a large suppression as observed in proton-nucleus (p-A) collisions at the SPS [7], the excitation function becomes rather flat. Already at that time, this demonstrated the importance of p-A collisions as a reference for the effects in AA systems.
However, indications for a much stronger suppression of the ψ(2S) have been observed in d-Au collisions at RHIC and more precisely established at the LHC, especially at backward rapidity (the nucleus-going direction) where the light-hadron multiplicity is the highest. These observations have been explained by final-state absorption on comovers [27], or, closely relate, by dissociation reactions in the QGP and hadronic phase of a fireball formed in these reactions [29,30] (both comover and thermal-suppression approaches feature comparable dissociation cross sections as well as fireball energy densities and timescales).
In the present paper we expand on our earlier results for d-Au collisions [30], by extending the kinetic rate-equation framework to p-Pb collisions at the LHC, including the calculation of p T spectra and rapidity dependencies. We construct an anisotropically expanding fireball based on initial asymmetries taken from Glauber model estimates of initial-shape fluctuations [31], which also allows us to compute charmonium elliptic flow. We recall that the rate equation approach necessarily includes regeneration contributions, which occur even in the presence of a single cc pair (sometimes referred to as "diagonal" regeneration or canonical limit). Their significance in p-Pb collisions has been suggested, e.g., in Ref. [28].
Our paper is organized as follows. In Sec. 2, we summarize the main components of the kinetic rate-equation/transport model developed for AA collisions and describe its extension to p-A collisions, in particular the anisotropic fireball evolution. In Sec. 3, we discuss our theoretical results for the nuclear modification factor as a function of centrality and p T , by first revisiting the d-Au system at RHIC, followed by 5.02 TeV and 8.16 TeV p-Pb collisions at the LHC, and compare to available experimental data from PHENIX, ALICE and LHCb. In Sec. 4, we discuss the v 2 results from our model in comparison to 8.16 TeV ALICE and CMS data. In Sec. 5, we summarize and conclude.
Our definition of backward (forward) rapidity in a p/d-A collisions follows the experimental convention of referring to the nucleus-going (proton-going) direction.

Transport Approach to Proton-Nucleus Collisions
The kinetic-rate equation approach developed in Refs. [32,33] is based on a space-momentum integrated Boltzmann equation, where N Ψ (τ ) denotes the time-dependent number of charmonium states Ψ = J/ψ, ψ(2S), χ c (1P). The transport parameters are the inelastic reaction rate Γ Ψ and the equilibrium limit (controlling the regeneration contribution), at a given temperature, T , of the ambient medium; V FB = V FB (τ ) denotes the timedependent fireball volume which we parameterize in Eq.(2.7) below, and f eq Ψ the thermal Bose distribution of the Ψ-state with degeneracy d Ψ . The charm-quark fugacity, γ c , is determined by charm-quark conservation in the fireball, where n op(hid) = n op(hid) (T ) denotes the density of all open (hidden) charm states in the system, and N cc is the total number of charm-quark pairs in the fireball which is calculated from the cc production cross section in pp (σ cc = dσcc dy ∆y with ∆y ≃ 1.8 for one fireball) at given collision energy N cc = σcc σ inel N coll , where N cc is the number of binary collisions at a given centrality of a p-A (or d-A) collision, and σ inel is the total inelastic cross section. The cc pair correlation volume, V co = 4 3 (R cc + v cc τ ) 3 (with initial cc pair radius R cc =1.2 fm and cc relative velocity v cc =0.6), represents the sub-volume that a pair can occupy after its production (sub-volumes of multiple pairs are merged once they overlap). Regeneration is only active for temperatures below the respective dissociation temperatures of each state Ψ (which are T c , 1.3 T c and 2 T c for the ψ(2S), χ c and J/ψ, respectively [33]). To account for incomplete charm-quark thermalization [34], the equilibrium limit has been corrected by a thermal relaxation time factor, where τ c is the thermal relaxation time of charm quarks, assumed to be constant [35]) at τ c =4 fm/c. Langevin simulations of charm quarks in p-Pb collisions [36] have found that a significantly larger relaxation time, by a factor of 3-5, is necessary to be compatible with the observed R pA of D-mesons; we therefore increase our previously employed relaxation time from 4 to 15 fm. The reaction rate Γ Ψ accounts for quasifree dissociation in the QGP and an extended set of dissociation reactions in the hadronic phase as constructed in Ref. [30] based on SU (4) meson exchange [37], supplemented by phase-space considerations for higher resonance states. Similar to Ref. [30], we increase the QGP dissociation rate of the loosely bound (or even unbound) ψ(2S) by a factor of 3 to account for non-perturbative effects on heavy-quark interactions in the QGP at moderate temperatures [38]. Uncertainties of the hadronic ψ(2S) dissociation rate are accounted for by increasing the baseline rate by up to a factor of 2. However, the variation of the hadronic rate has very little impact on the final R pA , due to a near compensation of dissociation and regeneration contributions. To schematically account for the effects of quantum evolution in the early stages of the charmonium evolution, we utilize formation times, τ form Ψ , for the different states (J/ψ, ψ(2S), χ c (1P )) that are assumed to have a range of 1-2 fm to reflect uncertainties associated with their binding energies. Their effect is rather small in semi-/central AA collisions, but becomes augmented in small systems due to shorter fireball lifetimes. The formation times not only influence (suppress) the magnitude of the early charmonium dissociation but also modify its p T dependence due to Lorentz time dilation implemented viã The time evolution of the fireball volume and temperature is constructed through an isotropic expansion with conserved total entropy, matched to the experimentally measured charged-hadron multiplicity in a given rapidity region of a nuclear collision system. The entropy density, s(T ), is calculated for a quasiparticle QGP and hadron resonance gas with a mixed phase at T c =180 MeV (in recent work [39] we have found that the use of a more modern lQCD-based EoS affects the temperature dependence of the fireball cooling (and thus the bottomonium kinetics) rather little; to keep consistency with our published charmonium results, we defer the EoS update to a future work).
Here, we extend our previously used cylindrical volume expansion to allow for a elliptic deformation in the transverse plane, where z 0 is the initial longitudinal size related to the formation time via τ 0 = z 0 /∆y, which we assume to be 0.8 fm, somewhat larger than in AA collisions to account for the reduced overlap density in the smaller p-A systems. The transverse radii in x-and y-direction are parameterized as The initial radii R 0 x and R 0 y are estimated from the eccentricity of the initial distribution of a Monte-Carlo Glauber event generator [31], with an initial transverse area A pPb ⊥ =πR 0 x R 0 y =7.8 fm 2 [31]. The surface velocities of the fireball are computed with a relativistic acceleration ansatz The parameters a x =0.34/fm and a y =0.13/fm are fixed in order to describe the lighthadron (pion, kaon, proton) p T spectra and v 2 including their mass splitting via the anisotropic blastwave formula, Eq. of a T ∼0.24/fm, relative to our default value of 0.1/fm in AA collisions, reflects the larger pressure gradients in p-A collisions [40]. We have checked that our total R pA results are rather insensitive to this value over a range of accelerations, a T =0.1-0.4/fm. Larger accelerations slightly reduce both the suppression and regeneration (compensating each other) due to shorter fireball lifetimes (most of the hot-matter effects happen at relatively early times where longitudinal expansion dominates). Blastwave fits of light-hadron spectra [41] extract average transverse velocities of up to ∼0.5 which would indeed require an transverse acceleration closer to ∼0.4/fm in our fireball framework. However, the blastwave fits in p-Pb collisions might be more sensitive, relative to AA collisions, to primordial hard components leaking into the fit range; thus a smaller transverse acceleration might be preferred. We therefore work with the values specified above unless otherwise stated. We also note that with the current fireball parameterization, the radii R x and R y cross at τ ∼2 fm, implying a transition to an in-plane deformation, cf. Fig. 1. We have checked using different ansätze that this is a robust feature dictated by a rapid build-up of the v 2 while approximately recovering light-hadron v 2 data in p-A collisions. After the crossing, the in-plane acceleration should become smaller than out-of-plane. While this feature is not explicitly guaranteed by our parameterization, the chosen parameter values lead to a sign flip of the anisotropic component of the acceleration, ∆a(τ ) = d∆v(τ ) , close to the crossing point of the radii. In any case, the net change in v 2 after this point is small and has very little bearing on our results.
Much like for different centralities in AA collisions, one can expect significant variations in the kinetic-freezeout temperature as a function of multiplicity in small systems: for a smaller total entropy the criterion that the mean-free-path is comparable to the fireball size (or inverse expansion rate) is reached at a larger particle density (or temperature). Guided by Ref. [41] we implement this effect by a centrality-dependent freezeout temperature as  The temperature evolutions following from this construction are summarized in Fig. 2 for 5.02 TeV and 8.16 TeV p-Pb collisions for different "centralities" (or rather, N ch ) at forward and backward rapidities.
To compute p T spectra within our approach, we decompose the rate equation into a solution for the primordial suppressed component, N prim , and the regeneration component, N reg . For the former, we solve a Boltzmann equation, where v Ψ denotes the charmonium velocity in the lab frame. The 2-dimensional vectors p T (p T , θ p ) and x T (r, θ r ) encode anisotropies in the transverse-momentum and coordinate plane, respectively, originating from different path lengths when traversing the elliptic fireball. The dissociation rate Γ Ψ ( p T , T (τ )) is evaluated in the medium rest frame where temperature is defined. In principle, there is a difference between the proper time τ in the rate equation (Eq.(2.1)) and the longitudinal proper time τ = √ t 2 − z 2 in the Boltzmann equation (Eq.(2.14)). However, since we place a single fireball at the appropriate rapidity for a given experiment, we do not correct its thermal rapidity width for time dilation. We also neglect the Lorentz contraction effects in the transverse volume expansion when solving the Boltzmann equation. The maximal effect at the surface with v T ∼0.5 is less than 15%, which is well within the uncertainty range of our acceleration parameter of 0.1-0.4 /fm.
For the regeneration yield in the elliptic fireball, we use the yield obtained from the rate equation and approximate its p T -spectrum by an anisotropic blastwave description [42,43], with α= p T T sinh ρ(r, θ r , τ ), β= m T T cosh ρ(r, θ r , τ ) and transverse mass m T = p 2 T + m 2 Ψ . The transverse-flow rapidity, ρ(r, θ r , τ )=tanh −1 (v ⊥ (r, θ r , τ )), is evaluated in terms of the fireball expansion velocity profile, v ⊥ (r, θ r , τ ) = (r/R max (θ r , τ ))v s (θ r , τ ), with surface velocity Since the surface velocity on the semi-minor axis (x-direction, with θ b =θ r =0) is larger than on the semi-major axis (y-direction, θ b =θ r =π/2), it generates a positive v 2 . The surface radius depends on the coordinate angle θ r and represents the boundary of the fireball, while characterizes the direction of the medium flow perpendicular to the fireball boundary.
To compute the denominator of the p T -dependent nuclear modification factor, and as an initial condition to the Boltzmann equation, we need the initial charmonium phase space distributions. We assume a factorization into transverse-momentum and coordinate space. For the p T distribution we employ an ALICE parametrization [44,45] of the spectra in pp collisions of the form with A=3.73(3.70), B=3.81(5.10) for J/ψ (ψ(2S)). The initial coordinate distribution is assumed to be a Gaussian, so that on the initial fireball boundary f = f 0 /e, and the enclosed elliptic area is equal to the initial transverse area (πR 0 x R 0 y ) within which all initial cc pairs are assumed to be produced (controlled by the normalization f 0 ).
The CNM effects are implemented in two steps. We first estimate the magnitude of the reduction (or enhancement) of the cc and Ψ yields from (anti-) shadowing using the EPS09-LO and EPS-NLO framework [46][47][48] at given rapidity and collision energy (defining an error band encoded in our final results). We then approximate the p T dependence of the CNM effects by a Gaussian broadening to represent both the original Cronin effect as well as the p T -dependence of shadowing, where a gN is the broadening per unit path length and L(b) the mean path length of the gluon in p/d-A collisions before fusing into a Ψ state [49]. It can be calculated as: for a nucleus A and proton/deuteron B, r A = | r − b 2 | and r B = | r + b 2 |. In the above expres- . (2.24) In the limit of zero absorption, it can be further simplified as which is used for the evaluation of the p T broadening. We associate the CNM effects for the p T dependence at the LHC with an effective broadening parameter of a gN =0.1-0.2 GeV 2 /fm reflecting the EPS09-LO vs. NLO uncertainty at backward rapidity, and a gN =0.2-0.4 GeV 2 /fm to represent the steeper trend and uncertainty from CGC calculations at forward rapidity [25,50]. At mid-rapidity, we take an intermediate range of a gN =0.1-0.3 GeV 2 /fm. The elliptic fireball allows the investigation of momentum anisotropies from final-state interactions. After obtaining the anisotropic spectra, dN AA /d 2 p T (p T , θ p ), from the pri-mordial and regeneration components, the elliptic flow coefficient is readily calculated as (2.26)

Nuclear Modification Factors for J/ψ and ψ(2S)
We are now in position to calculate the nuclear modification factors for charmonia in d-Au(0.2 TeV) collisions at RHIC (Sec. 3.1) and in p-Pb (5.02,8.16 TeV) collisions at the LHC (Sec. 3.2). The cross section inputs will be specified in the respective sections.

Deuteron-Gold Collisions at RHIC
Compared to our previous studies of d-Au collisions at RHIC [30], we here implement the updates as described in the previous section to ensure consistency with the new developments for p-Pb collisions. In particular, the fireball is extended to elliptic geometry, and has a smaller initial transverse area guided by the updates for p-Pb at the LHC described above; with a deuteron size approximately twice the proton size, and an inelastic N N cross section at RHIC of 2/3 of that at the LHC, we have A dAu As a result, the initial temperature in central d-Au now reaches T 0 ≃ 245 MeV. While this increases the hot-matter suppression, it slightly enhances the escape effect counter-acting the former. We also include regeneration contributions (neglected in Ref. [30]) which contribute up to ∼0.05 at the R dA level and also counter-act the increased hot-matter suppression.
The input cross sections remain unchanged, with dσ J/ψ dy =0.75 µb [51] for the J/ψ and dσcc dy =123 µb [52,53] for all cc pairs. Cold-nuclear-matter effects are associated with EPS09 LO parton shadowing [27,46] for both charmonium and open-charm production, whose centrality dependence we mimic by using a nuclear absorption cross section of σ abs =0-2.4 mb, as before [30]. The model calculations are compared with PHENIX data [8] in Fig. 3. Fair agreement with experiment is found, very similar to our previous results [30].

Proton-Lead Collisions at the LHC
In addition to the information specified in Sec. 2, we here quote the input cross sections as determined from pp data at the LHC. For the J/ψ we use dσ dy =3.0 and 3.6 µb at backward (-4.46<y<2.96) and forward (2.03<y<3.53) rapidity, respectively, at 5.02 TeV [54,55], and dσ dy =3.9(4.7) µb at backward (forward) rapidity at 8.16 TeV [21], and for cc pairs dσ dy =0.51(0.61) mb at backward (forward) rapidity at 5.02 TeV and dσ dy =0.66(0.80) mb at 8.16 TeV. This amounts to a fixed J/ψ over cc ratio of 0.58 % (as in our previous work [30,56]). The charged-particle multiplicity determining the total entropy of the fireball in the respective rapidity regions is extracted from Refs. [57,58] at 5.02 TeV and guided by Ref. [59] for 8.16 TeV. In the following, we first discuss the centrality dependence (Sec. 3.2.1) and then the transverse-momentum dependence (Sec. 3.2.2) of J/ψ and ψ(2S) production in 5.02 and 8.16 TeV p-Pb collisions.   . Centrality-dependent R pA for J/ψ (red bands) and ψ(2S) (blue bands) in 5.02 TeV p-Pb collisions, compared with experimental data [14][15][16]. The left (right) panel is for backward (forward) rapidity. The bands are due to (anti-) shadowing from EPS09 LO/NLO [47,48] at forward (backward) rapidity, as illustrated by the black bands which do not include final-state effects.

Centrality Dependence
In determining the centrality of a p-Pb collisions, we adopt the charged-particle multiplicity and relate it to the average binary collision number following Refs. [14,[57][58][59]. In Fig. 4, our J/ψ and ψ(2S) calculations are compared to 5.02 TeV ALICE data. The black bands show only the CNM effects, bounded by the anti-/shadowing obtained from EPS09-LO and EPS09-NLO calculations [47,48] for both charmonia and open charm; as for the RHIC case, the centrality dependence of shadowing is mimicked by a nuclear absorptiontype behavior, while for anti-shadowing we employ a parameterization of the pertinent lines shown in Fig. 3 of Ref. [27]. The CNM effects dominate the uncertainty bands at forward rapidity (charmonium formation time effects contribute ∼25%); the uncertainty bands at backward rapidity are entirely due to formation time effects (the same applies to Fig. 5). The shadowing-only bands already describe the J/ψ data quite well. A moderate hot-matter suppression of the J/ψ, together with a small regeneration contribution of about 0.05 (in units of the R pA ), generate additional suppression which leads to a slight underestimation of the backward-rapidity data but is compatible with the forward-rapidity data. For the ψ(2S) the much larger suppression in the hot fireball is, however, essential to approximately describe the suppression observed at both forward and backward rapidity. In Fig. 5, we compare our J/ψ and ψ(2S) calculations to 8.16 TeV ALICE data. There is a similar but slightly larger suppression compared to 5.02 TeV for both J/ψ and ψ(2S). We see quite some discrepancy with the data for peripheral collisions at backward rapidity, but fair agreement with the data at forward rapidity.

Transverse-Momentum Dependence
Our results for the p T dependence of charmonia, calculated as described in Sec. 2, are summarized in Figs. 6 and 7 for minimum-bias (MB) p-Pb collisions at 5.02 and 8.16 TeV, respectively. We recall that an additional uncertainty arises through the p T dependence of shadowing, which is incorporated into the theoretical bands conservatively as the maximum uncertainty from all effects. For the J/ψ, the calculated R pA (p T ) at backward rapidities at both 5.02 and 8.16 TeV exhibits a slight depletion at low p T followed by a mild maximum structure around p T ≃5-6 GeV, largely caused by the nuclear p T broadening. These trends become more pronounced at forward rapidity due to the generally increased strength of the CNM effects. Overall, the calculations are in agreement with ALICE and LHCb data within the theoretical and experimental uncertainties at both energies. The predictions for the ψ(2S) reflect the stronger suppression already observed in the centrality dependence. Relative to the J/ψ, most of the extra suppression is in the low-p T region where the hot-matter effects are most pronounced while at higher p T , formation time effects mitigate the suppression.

J/ψ and ψ(2S) Elliptic Flow
We finally turn to our calculation of the charmonium elliptic flow at 8.16 TeV, where data have recently become available [22,23]. The primordial component acquires a positive v 2 from the path length differences of the charmonium traversing the elliptic fireball, while the regeneration component acquires its v 2 from the anisotropic flow in a blastwave description. The primordial v 2 typically acquires values of 1-2%, while the v 2 of the regeneration component is much larger. However, since the latter, as mentioned above, is limited to R pA contributions of around 0.05-0.10, its weight in the total v 2 is small. Our results shown in Fig. 8 predict a small v 2 of up to ∼2% for the J/ψ, and a larger value of up to ∼5% for the ψ(2S), in high-multiplicity (most central) p-A collisions. We have checked tested that the maximal J/ψ v 2 generated from different versions of the fireball parametrization does not exceed 2%, essentially limited by the constraints from the initial eccentricity and the light-hadron v 2 . The near-zero result for the predominantly primordial component of the J/ψ is a direct consequence of its small hot-matter suppression (and regeneration): if it does not interact significantly, it cannot sense the spatial (or momentum) anisotropies in the fireball. This is also the reason why the v 2 of the ψ(2S) is much larger, since the hot medium effects on it are much larger. Since our J/ψ results clearly underestimate  Figure 8. Transverse-momentum dependent v 2 for J/ψ (red band) and ψ(2S) (blue band) at midrapidity in high-multiplicity p-Pb(8.16 TeV) collisions within the elliptic fireball model, compared to ALICE and CMS data [22,23].
the experimental data, we must conclude that the the observed v 2 cannot originate from final-state interactions alone. The similar v 2 at backward and forward rapidities (which have rather different multiplicities) is also in line with this conclusion. One last caveat we can think of are elastic interactions of the J/ψ (and ψ(2S)) in the expanding medium, which we have not accounted for. Very little is known about such interactions, and, in principle, one does not expect them to be large due to the parametrically smaller size of the J/ψ compared to light hadrons, while for the ψ(2S), due to its small binding, almost any interaction can lead to break-up.

Conclusion
In the present work, we have extended our transport approach for in-medium quarkonia in heavy-ion collisions to calculate J/ψ and ψ(2S) production in small collision systems at RHIC (d-Au) and the LHC (p-Pb). Cold-nuclear-matter effects estimated from nuclear parton distribution functions are combined with final-state effects treated within a rate-equation framework for an expanding fireball including dissociation and regeneration reactions in the QGP and hadronic phase. Our calculations provide a generally fair description of the measured centrality and transverse-momentum dependent nuclear modification factors measured in different rapidity regions, which differ in their CNM and hot-nuclear matter effects (some tension with data was found in the 8.16 TeV backwardrapidity R pA (N coll )). This supports an interpretation where the J/ψ observables are mostly dominated by CNM effects while the loosely bound ψ(2S) is subject to substantial suppression in the hot fireballs with initial temperatures of about 200-300 MeV and lifetimes of up to 4 fm. We also investigated the elliptic flow of J/ψ and ψ(2S). In our setup, a nonzero v 2 results entirely from final-state interactions in the elliptic fireball. Since the final-state suppression (and regeneration) especially for the J/ψ is small, which is compatible with the small hot-matter effects on the R pA , the resulting v 2 is also small, not more than 2% (and larger, up to 5%, for the ψ(2S)); this disagrees with the large signal observed in the LHC data. We are therefore forced to conclude that this signal must be in large part due to initial-state (or pre-equilibrium) effects not included in our approach. This situation appears to be part of a bigger picture where the nuclear modification factor of hadrons, e.g., D-mesons, shows little deviation from one while the v 2 is appreciable.