A simple model to include initial-state and hot-medium effects in the computation of quarkonium nuclear modification factor

Quarkonium suppression is one of the more useful observables to obtain information about the hot medium created in ultrarelativistic heavy-ion collisions. In this manuscript, we discuss a simple way to implement both the initial-state effects and the hot-medium evolution, and to compute the quarkonium nuclear modification factor if the survival probability for a bound state at a given energy density is known. Based on the Glauber model, the temperature of the evolving medium and the centrality dependence of the nuclear modification factor will be analysed.


I. INTRODUCTION
Quarkonium suppression was proposed in [1] as a signal of a formation of a quark-gluon plasma. Over the years, the models for the description of quarkonium evolution in a hot medium have been refined [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16] and our theoretical description is now very different to the original idea. However, to obtain a realistic description of quarkonium in a heavy-ion collision, it is not enough to understand how quarkonium interacts with a medium. In the following, we mention just some of the relevant issues: • We need to take into account the time evolution of the medium itself. Relativistic hydrodynamics is one of the most common frameworks [17] as it gives an excellent description of the observables related to the bulk of the medium. If we assume that the nucleus properties are homogeneous in the transverse direction, we obtain the Bjorken evolution [18]. In any case, we need to add as input some information about the initial state of the medium, usually the initial energy density at any point in the transverse plane, and the initial time at which a hydrodynamic description of the medium is valid.
• Naively, we would expect that we could predict what would happen to quarkonium in a heavy-ion collision if the medium would not modify quarkonium's population from extrapolating proton-proton data. However, this is very far from reality. It is quite non-trivial even to predict proton-nucleus collisions, in particular due to the presence of initial state effects that lead to the modification of the nuclear Parton Distribution Functions (nPDFs). Several models exist for this task [19][20][21][22].
Here we propose to use the model of [19], based on Pomeron interactions. In this model, the origin of the modification of the nuclear ratios is related to multiple scattering with nucleons.
This results in a reduction of the corresponding cross sections, commonly called shadowing. This simple model offers the advantage of being completely analytical and provides a reasonable description of proton-nucleus data.
Our aim in this manuscript is to provide a simple method to compute the nuclear modification factor in heavy-ion collisions, taking into account initial nuclear matter and hot medium effects. To do so, we will assume that the following information is available: • The initial state of a heavy quark pair after a proton-proton collision.
• The probability that this heavy quark pair forms a given bound state provided that it interacts with an evolving medium which initially has an energy density ǫ i .
First, we will compute how the initial temperature depends on the position on the transverse plane. This temperature depends on the medium's energy density, that we assume proportional to the particle multiplicity, i.e. to the number of pions. We describe the multiplicity using the model of [19] which takes into account nuclear shadowing effects. We argue that to a very good approximation, this is equivalent to assuming that the initial energy density at a given point in the transverse plane is proportional to the density of participant nucleons at that point. This statement is true at least at nowadays collinear energies. To illustrate this, we will show our results for the temperature behaviour in 3 different cases, i.e. considering it to be proportional to the number of participants, to the number of collisions and to the number of pions, once shadowing effects are taken into account. Then, we will illustrate how to compute the quarkonium's nuclear modification factor in heavy-ion collisions including initial and hot nuclear matter effects. We remark that our aim is not to model quarkonium interaction with a medium but to provide a prescription to include initial effects and temperature evolution for a quarkonium computation in a hot environment.
This prescription has the advantage that the effects of shadowing can be incorporated in a straight-forward way using analytical formulas.
The outline of the manuscript is as follows. In section II, we review some aspects of the Glauber model useful to set the scope and notation of our computation. In section III, we review the initial state model of [19]. In section IV, we study the initial energy density and the temperature of the medium. In section V, we show how to implement our model for the computation of the nuclear modification factor. In section VI, we illustrate our prescription with some examples in which we include both hot and initial nuclear matter effects on quarkonium. Finally, in section V, we give our conclusions.

II. GLAUBER MODEL REVIEW
For completeness, we review in this section the Glauber model [23] in the optical approximation [24]. This is needed in order to fix the notation. Many reviews of the Glauber model can be found in the literature, in our case we follow closely the description given in [25].
The optical Glauber model is used to extrapolate the known properties of the collisions between nucleons to describe collisions of large nuclei made of many nucleons. For all practical purposes within this manuscript, we take the following physical assumptions: • The eikonal approximation is used for the nucleons inside a nucleus. This means that for the purposes of studying the nucleus-nucleus collision we can consider that the nucleons do not have transverse motion. This approximation is correct up to subleading corrections of order Λ QCD √ s .
• The interaction between nucleons is local, therefore, two nucleons of different nucleus can only interact if they are at the same point in the transverse plane during the collision.
• Collisions between nucleons are independent processes. In other words, the probability of collision between two nucleons does not depend on whether the rest of the nucleons have collided or not.
The function that contains all the required information of the nucleus is the density of nucleons inside a nucleus, which is assumed to be proportional to the nuclear charge density.
We focus on the case of a spherically symmetric nucleus and we use a Woods-Saxon profile where for the case of lead c = 6.62 fm and ξ = 0.546 fm (These values are taken from [26]).
A is the nucleon number and ρ 0 is chosen such that An analytical expression can be found for ρ 0 where Li is a polylogarithm. From the density of nucleons we can compute the thickness function, which informs of the density of nucleons at a given position in the transverse plane for any value of the longitudinal coordinate Since we are interested in the collision of a nucleus with nucleon number A with another with nucleon number B, it is useful to define the overlap function, i.e. the density of pairs of nucleons from different nuclei but at the same point in the transverse plane -so that they might interact. This can be defined for a given value of the impact parameter b as Note that in the previous equation we use a coordinate system for the transverse plane such that s = 0 corresponds to the center of the colliding region when the two nuclei are identical.
This is not the standard choice in the literature, where it is common to choose an origin of coordinates which coincides with the center of one of the nuclei. Our motivation is that we expect many of the properties of the collision to exhibit a central plateau and we want a coordinate system whose origin is closed to the center of the plateau.
Until now we have only described the geometry of the collision. Another required input is the cross section for inelastic collision between two nucleons, which in order to simplify the notation we represent with σ. At the LHC, the value σ = 64 mb is taken for √ s = 2.76 TeV collisions while the value σ = 70 mb is taken for √ s = 5.02 TeV collisions. With this we can compute the probability of having n nucleon collisions. Since we assume that the collisions are independent processes it must follow a binomial distribution. We consider the number of n pairs that can be made out of A nucleons from one nucleus and B from the other, then we consider the probability that these n pairs collide while the rest of possible pairs do not.
This reasoning leads to the following expression Then the probability that there is at least an inelastic collision is and the mean value of the number of collisions is Another quantity which is used to characterise heavy-ion collisions is the number of participants, i.e. the total number of nucleons that collide at least once. To compute it we can use the fact that, in the optical Glauber model, collisions are assumed to be independent.
Thus we can compute the number of participants in the following way: • Take one nucleon from the nucleus with mass number A.
• Compute the probability that this nucleon has an inelastic collision with the other nucleus. To do so we can use eq. (7) adapted to the pB case (P (s) inel,pB = 1 − ).
• Then we sum over the distribution of nucleons in the nucleus with mass number A.
• We repeat the process but now studying the probability that a nucleon of nucleus B interacts with nucleus A.
Therefore, we arrive at the following expression In many cases, experimental results are given in terms of centrality classes. It is assumed that the multiplicity increases as the impact parameter decreases. Then, for example, the 10 % of collisions with higher multiplicity corresponds to the 10 % more central collisions.
Therefore, the 0 − 10 centrality window includes collisions from b 0 = 0 up to another impact parameter that we call b 10 . The numerical value of b 10 can be computed in the following way. Let us assume that the multiplicity is proportional to the inelastic cross section, given in eq. (7). It is convenient to simplify eq. (7) using that AB ≫ 1 It follows that the total inelastic cross section is Then b 10 can be determined solving the following equation Analogously, the centrality window x−y includes collisions in the range of impact parameters between b x and b y such that

III. INITIAL-EFFECT MODEL
As it is commonly assumed, the quarkonium production, primarily formed by gluon-gluon fusion, will suffer the shadowing effects that modified the parton distribution functions in a nucleus and affect global multiplicities. Several models, mostly based on data parametrization at a given Q 2 followed by Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution, exist [20,21]. For our purpose, it is more useful to implement the analytical shadowing model based on Glauber-Gribov theory [19]. In this model, the interaction of two high- In order to relate diffraction on nucleons with nuclear shadowing, the γ * nucleus cross section can be expanded in a multiple scattering series containing the contribution from γ * A is just equal to Aσ γ * nucleon , the Glauber elastic contribution. The amount of the shadowing corrections corresponding to two-scattering contribution obeys Eq. (14) represents the first contribution to nuclear shadowing originated by double scattering with two target nucleons. Convoluted with the normalised nuclear profile function, it gives a negative contribution to the total cross section, σ . The ratio of the differential cross section for diffractive dissociation of the virtual photon over the total γ * nucleon cross section 1 σ γ * nucleon dσ D γ * nucleon dM 2 dt | t=0 that appears in the above equation is directly related to the ratio of the triple-Pomeron cross section over the single Pomeron and can be interpreted as the reduction due to the interaction among the partons.
Thus, in terms of triple-Pomeron contributions, eq. (14) reads where we have introduced the standard triple-Pomeron formula which involves the Pomeron intercept and its residue, 4π 1 , being g P pp the Pomeron-proton coupling and r P P P the triple-Pomeron coupling, both evaluated at t = 0, C = 0.023 fm 2 . The value of ∆ is related to the value of the Pomeron intercept α P (0) = 1 + ∆ = 1.13.
We consider that all intermediate states have the same structure and therefore we use the Schwimmer scheme to include higher-order rescatterings [28]. The suppression from shadowing in AB collisions for particles produced at rapidity y is obtained replacing the nuclear profile function T AB (b) by the integral on s of where is the triple Pomeron graph contribution, with the rapidity of the triple Pomeron vertex integrated up to y (where the trigger particle is produced), i.e. up to y max = y + ln ( √ s m T ). y is the center of mass rapidity of the produced particle, y > 0 for the projectile hemisphere (forward rapidity) and y < 0 for the target one (backward rapidity). This integration limit is related to the parton model for hard processes: for projectile A (target B), x A(B) = m T √ s e ±y . m T is the transverse mass of the particle, m T = m 2 + p 2 T . We take y min = ln ( Here m N is the nucleon mass and R A = 0.82 A 1/3 + 0.58 fm is the Gaussian nuclear radius. T A and T B are the nuclear profile functions for which, as mentioned above, we use a Woods-Saxon parametrization. Note that the amount of shadowing depends on the nature and the transverse momentum of the produced particle through its transverse mass. dy for the yields of produced particles in the central rapidity region at high energies. Comparison with data shows that this relation significantly overestimates the hadron multiplicity at LHC energies. The corrections due to shadowing effects, arising by the contribution of the triple Pomeron graphs, cure this discrepancy and change the dependence of the multiplicity which behaves that way proportional to A δ . The value of delta is a function of energy and it is equal to δ ≃ 1.1 at LHC energies. Within this approach, a good description of the multiplicities from the Relativistic Heavy-Ion Collider (RHIC) to LHC energies is obtained [19].

IV. ENERGY DENSITY AND TEMPERATURE OF THE MEDIUM
In order to get the energy density and the temperature of the medium, we will carry out the following approach. The energy deposition in the transverse plane corresponds to the transverse mass of the produced particles, so it is proportional to the multiplicities. As a first approach, in order to get the initial energy density one could take it as proportional to the density of wounded nucleons N part . While at low energies this is a common assumption, the second contribution is expected at RHIC and LHC energies, proportional to the number of collisions. However, at RHIC energies, the PHOBOS collaboration has checked at different energies (20, 60 and 200 GeV) that the total charged particle multiplicity exhibits an essentially linear dependence on N part [29]. Comparison of the averaged and scaled 200 GeV AuAu data with LHC data on multiplicities versus centrality shows a remarkable similarity in the shape of both distributions, at least up to peripheral collisions.
The fact that the shape of the normalised multiplicity distribution varies little with energy and stays almost constant up to TeV energies -even if the amount of hard processes, which scale with the number of binary collisions N coll , is expected to contribute significantly to particle production at LHC and should lead to a steeper centrality dependence-can be explained by a strong impact-parameter dependent shadowing of the nPDFs that limits this rise with centrality. This is illustrated in Fig. 1, where the initial temperature at the center of the colliding region, taken as proportional to the energy density, T 0 (s, b) ∝ (ǫ i (s, b)) show that, up to peripheral collisions, the simplification that consists on taking the energy density proportional to the number of participants is correct and offers the same result as the one based on the multiplicity. In the following, we will use the results based on the multiplicities, which includes shadowing, to determine the initial temperature. In practice, the numerical cost of determining the temperature considering pion shadowing is almost equal to that of assuming that the temperature goes like the local number of participants to 1/4.

V. APPLICATION TO QUARKONIUM SUPPRESSION
In this section, we discuss how the previous results could be used to compute the nuclear modification factor R AA of quarkonium in a heavy-ion collision. R AA is the quantity of a particular quarkonium state observed in a heavy-ion collision divided by the number that we would naively expect from proton-proton data. More precisely, in absence of any cold or hot nuclear matter effect any collision of nucleons would be approximately equivalent to a proton-proton collision. Then, since the number of collisions is proportional to T AB (b), R AB can be computed with where N AB HQ is the number of quarkonium states observed in heavy-ion collisions -for simplicity N HQ -and N pp HQ is the cross section of quarkonium in pp collisions. We can estimate N HQ as where S sh HQ corresponds to the density of binary collisions corrected by the amount of shadowing on quarkonium estimated according to the model discussed in section III and S med is the medium induced suppression. In this way, if we assume that S med = 1, we would be taking into account only initial nuclear matter effects. On the other hand, we can ignore initial nuclear matter effects by substituting S sh HQ (s, b) by T AB (s, b), i.e. taking F (y) in eq. (16) equal to 0. We name the quantity obtained by setting S med = 1, R I AA . Instead, if we substitute S sh HQ (s, b) by T AB (s, b) we obtain R T AA . Note that we are assuming that quarkonium is comoving with the medium. To consider finite velocity effects is out of the scope of this work and is left for future investigations. A quarkonium that is not comoving with the medium would be sensitive to the temperature of different medium regions at different times. The initial condition model we are using would still be valid in this situation, however, eq. (19) would not be valid anymore.
The effects of the hot medium are usually taken as S med (s, b) = f (T 0 (s, b)), where T 0 is the initial temperature of the medium. In other words, the medium suppression depends on the initial temperature at the point in which the quarkonium pair is created. Taking the initial energy density proportional to the number of pions created at a given point the calculation of the initial temperature is straightforward. A common assumption is where we have considered that the equation of state is that of a free gas. We assume that the energy density at a given point is proportional to the production of pions. Then, the temperature just after the collision is given by where, by construction, T 00 is the temperature in the center of the overlapping region for the most central collisions [30]. This is a quantity whose computation is out of the scope of this work. It can be adjusted using experimental results within a hydrodynamical framework [31,32]. The value for PbPb collisions at √ s = 5.02 TeV is estimated to be 502 MeV.
We are going to focus on computing the relation of the temperature at different points and values of the impact parameter with T 00 . We can compute the average initial temperature .
A good approximation to the initial temperature in a not too peripheral collision is given by T 0 (0, b). The reason is that T AB has a central plateau in which it is approximately constant.
We compare T 0 (0, b) and T HQ (b) in Fig. 2. We see indeed that to take the temperature at the center of the overlapping region as representative of the temperature of all the collision is a reasonable approximation, at least regarding the temperature seen by a quarkonium state.
In Table I  In what follows, we will show the applicability of our approach within several models of the quarkonium-medium interaction. In the previous sections we have discussed how R AA can be computed ignoring medium effects by taking eq. (19) in the case S med = 1. The computation of the survival probability of quarkonium in medium is a very active research topic in itself. In this section we discuss the results of two different scenarios. In the first one, we make use of a simplified model inspired by recent potential non-relativistic QCD computations [33,34]. In the second one, we use the gap model developed in [35] that implements the lattice static potential. . Note that, because Γ scales like T 3 , β is a constant that does not depend on b and s. Then we can rewrite, This is the formula that we are going to use in our computations of thermal effects in this subsection. α and β have been chosen such to qualitatively follow experimental results.
Let us remark again that our purpose here is not to make a state-of-the-art description of quarkonium in a medium but rather to have a simple analytic formula for quarkonium's survival probability that we can use to illustrate how to include cold nuclear matter effects in a thermal model within our approach.
In order to compute R AA we use eqs. (18), (19) and (21). In Fig. 3, we show our results for (which is considered to include our initial-state effects) and R T AA (see definitions in section V) using eq. (27) with α = 0.4 and β = 1.2. Although the model is very naive, the computation allows us to obtain some conclusions. Shadowing effects are significant, although smaller than medium effects. However, they have some significant qualitative differences. Shadowing effects have a milder dependence with the number of participants and they decrease slowly as this number goes to zero. It is also important to note that the total suppression is not the sum of the suppression due to shadowing plus the suppression due to in medium effects. This implies that the influence of shadowing on R AA depends on the specific model for quarkonium-medium interaction, at least if the aim is to go beyond qualitative considerations.
We remark again that the survival probability used in this subsection is a simplified model based on potential NRQCD (pNRQCD) arguments valid in the region 1 r ≫ T ≫ E. The parameters in eq. (27) were chosen in an ad hoc way to qualitatively reproduce experimental results on R AA once shadowing effects are taken into account. A more precise, rigorous and state-of-the-art application of the pNRQCD framework to the computation of R AA , ignoring shadowing effects, can be found in [34,39].

B. Gap model
In this subsection we aim to include shadowing effects in the model developed in [35].
More specifically, we focus on the scenario discussed in [35] in which lattice QCD data on the static potential is used to compute the decay width of Υ(1S). The emphasis of [35] was to highlight the importance of taking into account the finite energy gap between color singlet and color octets states. The decay width of quarkonium is suppressed compared to computations which do not take into account the finite value of the energy gap. This is a remarkable effect, especially at small temperatures. In [35], it was shown that suppression is overestimated by an approximate factor of two when the gap is ignored.
The S med that we are going to use in this section is where a = 22.9, b = 2170 MeV, t 0 = 0.6 fm, T f = 175 MeV and T 0 is the initial temperature.
Details about the derivation of this formula can be found in [35]. quantum system approach in the limit E ≫ Γ. The decay width takes into account the finite energy gap between singlets and octets. Lattice QCD data on the static potential is used to determine the wave function and the binding energies of quarkonium, both information is needed to compute the decay width. In [35], the decay width was computed for several temperatures. The results could be well fitted by the following function This formula was then used to derive eq. (28). Here we improve the previous computation by including shadowing effects. They influence R AA in two different ways. By modifying the value of T 0 in eq. (28) and by the shadowing effects included in S sh HQ at eq. (19). Our results are shown in Fig. 4. We note that R T AA is almost identical to the R AA results obtained in [35]. This is the expected result since the only difference is that in the computation of R T AA pion shadowing was taken into account while in the [35] it was assumed that the energy density scales like the local density of participants. This is, as we discussed, a very good approximation at LHC energies. Note also that R CN M AA comes from the same formulas in Figs. 3 and 4, the only difference is due to statistical fluctuations in the Monte Carlo integration done to obtain the results. Focusing on R AA , including both shadowing and medium effects, we observe that shadowing effects are practically the only source of suppression when the number of participants is below 50. Our interpretation is that at very small temperatures the energy gap suppresses the decay width very strongly, leaving shadowing as the only source of suppression. Finally, the results in Fig. 4 reinforce the observation that the overall suppression is not simply the sum of the suppression due to shadowing plus the suppression due to thermal medium effects.

VII. CONCLUSIONS
In this manuscript we have discussed how shadowing corrections can be easily included in the computation of quarkonium suppression. Our main assumptions are that the medium induced suppression can be described by a survival probability, the validity of the shadowing model discussed in [19], that the initial temperature scales as 1/4 of the initial energy density and that this quantity is dominated by the pions' contribution. We note that our method could in principle be used in combination with a Markovian quantum evolution of quarkonium in a medium (as that in the recent paper [39]) as long as recombination effects are ignored.
One of our main focuses is the temperature felt by a quarkonium state formed at a given point in the medium. We found that to assume that the energy density is proportional to the density of participants gives almost the same result as considering a more realistic model of pion shadowing, at least regarding the initial temperature. Another interesting observation is that the average temperature seen by a quarkonium state is approximately equal to the temperature in the center of the overlapping region.
Finally, we have applied our methods to two different models of the medium-quarkonium interaction. The first is a model inspired by recent pNRQCD results in which the medium sees quarkonium as a small color dipole and the temperature is much larger than the binding energy. The second is one of the scenarios considered in [35]. In both cases we found that shadowing effects are smaller than medium effects but of the same order of magnitude.
Shadowing has a milder dependence on the number of participants than medium induced suppression. Finally, it is interesting to observe that the overall suppression is not simply the sum of the suppression due to shadowing plus that due to medium effects.
In principle, our method could also be used with models in which medium effects are a mixture of quark-gluon plasma and hadron gas effects, as long as these effects can be encoded in a survival probability.