Aspects of relativistic heavy-ion collisions

The rapid thermalization of quarks and gluons in the initial stages of relativistic heavy-ion collisions is treated using analytic solutions of a nonlinear diffusion equation with schematic initial conditions, and for gluons with boundary conditions at the singularity. On a similarly short time scale of $t \le1$ fm/$c$, the stopping of baryons is accounted for through a QCD-inspired approach based on the parton distribution functions of valence quarks, and gluons. Charged-hadron production is considered phenomenologically using a linear relativistic diffusion model with two fragmentation sources, and a central gluonic source that rises with $\ln^3(s_{NN})$. The limiting-fragmentation conjecture that agrees with data at energies reached at the Relativistic Heavy Ion Collider (RHIC) is found to be consistent with Large Hadron Collider (LHC) data for Pb-Pb at $\sqrt{s_{NN}}= 2.76$ and $5.02$ TeV. Quarkonia are used as hard probes for the properties of the quark-gluon plasma (QGP) through a comparison of theoretical predictions with recent CMS, ALICE and LHCb data for Pb-Pb and p-Pb collisions.


I. INTRODUCTION
This article covers aspects of relativistic heavy-ion collisions in the energy regions reached at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC). Starting with the thermalization of partons in the initial stages, stopping of the incoming baryons is discussed, followed by charged-hadon production and the modification of quarkonia in the hot quark-gluon plasma (QGP), with an emphasis on bottomonia that provide clear signals of QGP formation. The overall approach is phenomenological and in close comparison with data. In some cases such as stopping and quarkonia, nonequilibrium-statistical considerations are merged with quantum chromodynamics (QCD). The work contains both new developments and re-views of our previously published work in a new context.
The fast local thermalization of partons in the initial stages of relativistic heavy-ion collisions is a sufficient condition to apply hydrodynamic descriptions [1] of the subsequent collective expansion and cooling of the hot fireball that is created in the collision. Typical local equilibration times for gluons are about 0.1 fm/c [2], with initial central temperatures of the order of 500-600 MeV reached in a Pb-Pb collision at √ s NN = 5.02 TeV at the Large Hadron Collider [3]. Thermalization times for quarks are typically by an order of magnitude larger [4] due to the smaller color factor and the different statistical properties (Pauli's principle). Whereas the thermal Bose-Einstein/Fermi-Dirac distributions and also the initial distribution for quarks are known, plausible assumptions are being made for the primordial gluon distribution before the onset of the collision.
Studies such as Ref. [5] have found that local thermalization may not be a necessary condition for the appli- * g.wolschin@thphys.uni-heidelberg.de cability of hydrodynamics to relativistic heavy-ion collisions. This would imply that hydrodynamics could be a valid approach away from local equilibrium [6] -but still, it remains important to investigate how and on which timescale local thermalization is achieved. Obviously, an investigation of the local thermalization of quarks and gluons makes sense only in a weakly coupled description based on an effective kinetic theory that relies on the Boltzmann equation. It is recognized that a strongly coupled paradigm built on the anti-de Sitterconformal field theory (AdS/CFT) correspondence [7] discovered in the investigation of D-branes in string theory [8] may be relevant at RHIC and LHC energy scales such that partons would not be the relevant degrees of freedom any more. In this work, however, I rather assume quarks and gluons to be well-defined and long-lived excitations in QCD at temperatures close to the critical value.
What is not known precisely is the development with time towards local equilibrium. An easy way to approximate the time evolution is given by the linear relaxationtime approximation (RTA), which provides a simple analytic solution of the problem by enforcing a linear transition from the initial to the local equilibrium state, but does not properly account for the known nonlinearity of the system. A more ambitious approach for bosons is to numerically solve gluon transport equations that include the effect of Bose statistics, as has been done in Ref. [9] and subsequent works. Initially, only elastic scattering was considered with the possibility of gluon condensate formation due to particle-number conservation in overpopulated systems, but it was recognized that inelastic, number-changing radiative processes cannot be neglected [10], and hinder the formation of a Bose condensate.
Whereas also other numerical calculations relying on a quantum Boltzmann collision term to account for the initial local equilibration are available such as Ref. [11], it is of interest to have an exactly solvable analytic model to better understand the physics of the fast equilibration. A corresponding nonlinear boson diffusion equation (NBDE) has been presented in Ref. [4] and solved for a simplified case that did, however, not yet consider the singularity at the chemical potential µ < 0. The nonlinear partial differential equation preserves the essential features of Bose-Einstein statistics that are contained in the collision term. In particular, the thermal equilibrium distribution emerges as a stationary solution and hence, the equation appears suitable to model the thermalization of gluons in relativistic collisions. It is used in this work to generate new exact solutions for the time-dependent gluon equilibration problem that include boundary conditions at the singularity. Regarding fermionic thermalization, the corresponding nonlinear fermion diffusion equation is easier to solve [4,12] because no singularity appears, and the analytic solutions will be reviewed.
On an equally short time scale as the local equilibration, the incoming baryons with energies available at RHIC or LHC are stopped: The system is slowing down, essentially through collisions of the incoming valence quarks with soft gluons in the respective other nucleus. Various models to account for this process and its energy dependence have been developed, in particular in Refs. [13,14] and related works, which are relying on the appropriate parton distribution functions (PDFs) and hence, on QCD, yielding good agreement with the available net-proton (proton minus antiproton) stopping data. Different from the nonequilibrium-statistical approach to initial thermalization, such models do not consider a time dependence. However, by using the rapidity distribution calculated from the PDFs and the initial valence-quark distribution, one can, in addition, account for the time development from the initial to the final distribution with an appropriate fluctuation-dissipation relation.
The relevant bulk properties of relativistic heavy-ion collisions mostly arise from charged-hadron production. In many of the available macroscopic and microscopic models, hadronization occurs from the fireball at the phase boundary between the QGP and the hadronic phase. However, it has been proposed [15,16] that the production of charged hadrons from the fragmentation sources at larger rapidity values is also relevant in the overall distributions and should be treated separately from the fireball source. At midrapidity, these contributions are relatively small, but become relevant more forward or backward. Corresponding new results of the phenomenological three-source relativistic diffusion model (RDM) are compared with recent LHC pseudorapidity data on charged-hadron production in 5.02 TeV Pb-Pb, which are also shown to be consistent with the limiting-fragmentation hypothesis that had been found to agree with the hadron production data at RHIC energies [17][18][19].
Regarding more direct evidence for transient QGP formation in relativistic heavy-ion collisions, jets may pro-vide the most direct manifestation for quarks and gluons in the system. In particular, the suppression of awayside jets in the hot QGP had already been predicted by Bjorken [20], and confirmed experimentally by the STAR collaboration [21]. Meanwhile, jet suppression at LHC energies has been investigated in detail (e. g. Ref. [22]), and it has been shown how strong final-state interactions cause high-p T jets to lose energy to the plasma.
In this work, however, I consider quarkonia as another indicator for the properties of the quark-gluon plasma such as its initial central temperature T i . Quarkonia are bound states of heavy quark-antiquark pairs that can be formed in initial hard partonic collisions. In the original prediction for a suppression of the J/ψ yields in the presence of a QGP, only the medium-effect on the real part of the quark-antiquark potential was considered [23]. It has later been realized that due to the presence of the hot medium, the potential has an imaginary part [24]. Optical potentials had also been used in the theory of nuclear reactions to account for channels that are not treated explicitly. In case of quarkonia, the imaginary part causes their dissociation, in addition to melting of the quarkonia states because the real potential is screened. It is possible to treat quarkonia dissociation by thermal gluons separately from the imaginary part [25].
In case of charmonium at LHC energies, statistical recombination of charm and anticharm quarks turns out to be important, but it is not possible to separate the process from the dissociation in the QGP. One may therefore concentrate on the heavier bottomonium system, where recombination is much less pronounced. Through detailed investigations of the transverse-momentumand centrality-dependent suppression of the spin-triplet Υ(1S, 2S) states in Pb-Pb collisions at LHC energies and comparisons with CMS data, we can deduce QGP properties such as the initial central temperature T i and study its dependence on the center-of-mass energy.
In asymmetric collisions such as p-Pb, the situation is quite different from symmetric systems, because most of the system remains cold due to the much smaller overlap. Cold nuclear matter (CNM) effects therefore provide a certain understanding of the measured quarkonia modifications, but a complete agreement with the available data remains impossible -unless one also considers the hot QGP zone that is still produced, even though it is initially considerably less extended as compared to symmetric systems. Indeed the bottomonia dissociation in the hot QGP provides a clue for the interpretation of the data in asymmetric collisions as well.
The paper follows the above-mentioned series of topics: In Section 2, the nonlinear diffusion equation for gluons and quarks is solved explicitly to account for the fast local thermalization at the beginning of relativistic heavy-ion collisions. Stopping is considered in Section 3, hadron production and limiting fragmentation in Section 4, bottomonia modification in the medium in Section 5. The conclusions are drawn in Section 6.

II. FAST THERMALIZATION OF GLUONS AND QUARKS: AN ANALYTIC NONLINEAR MODEL
For a given initial nonequilibrium gluon distribution at t = 0, solutions of the nonlinear boson diffusion equation describe the time-dependent equilibration towards the thermal distribution with the local temperature T . In Ref. [4] such solutions were calculated with the free Green's function. Whereas this accounts for local thermalization in the ultraviolet (UV) with the corresponding equilibration time τ eq , in the infrared (IR) the populations decrease due to diffusion into the negative-energy region. To avoid such an unphysical behaviour, one has to consider the boundary condition at the singularity |p| = p = ǫ = µ with the chemical potential µ < 0, and the corresponding bounded Green's function in the solution of the NBDE. With this Green's function, gluon populations indeed attain the Bose-Einstein limit also in the infrared for nonequilibrium initial conditions that include the singularity.
The nonlinear model and the solution of the combined initial-and boundary-value problem are first briefly reviewed. Subsequently, the thermalization problem is solved for a schematic initial gluon distribution that characterizes the relativistic collision at t = 0. Adding the boundary condition at the singularity, n(ǫ = µ < 0, t) → ∞, the time-dependent partition function that includes initial and boundary conditions is obtained using analytic expressions for both, the bound Green's function, and the function that contains an integral over the initial conditions. The resulting occupation-number distribution function n(ǫ, t) is calculated, and it is shown to approach the equilibrium distribution both in the UV and in the IR.
Before we proceed to the nonlinear model, it is useful to consider a linear time-dependent transition from the initial distribution to the thermal distribution with the chemical potential µ < 0 in a finite boson system in the relaxation-time approximation (RTA), ∂ n rel /∂t = (n eq − n rel )/τ eq , with solution Time-dependent RTA results for gluons are shown in Fig. 1 for t = 0.02, 0.08, 0.15, 0.3, and 0.6 fm/c. The thermal distribution with initial central temperature T = 513 MeV as inferred in Ref. [3] for central Pb-Pb collisions at √ s N N = 5.02 TeV is approached linearly, the discontinuities at ǫ = Q s persist.
For a more realistic account of thermalization, one needs to consider the inherent nonlinearity of the sys-  Starting from schematic initial conditions Eq. (13) in the cold system at t = 0 (box distribution with cut at ǫ = Qs = 1 GeV), a Bose-Einstein equilibrium distribution with temperature T ≃ 513 MeV (solid curve) is approached. Time-dependent single-particle occupation-number distribution functions are shown at t = 0.02, 0.08, 0.15, 0.3, and 0.6 fm/c (decreasing dash lenghts). The lower dotted curve is Boltzmann's distribution.
tem. I had derived a nonlinear partial differential equation for the single-particle occupation probability distributions n ≡ n th (ǫ, t) from the bosonic/fermionic Boltzmann collision term in Ref. [4]. The transport coefficients in this nonlinear boson diffusion equation (NBDE) depend on energy, time, and the second moment of the interaction. They incorporate the complicated many-body physics. The drift term v(ǫ, t) is responsible for dissipative effects, the diffusion term D(ǫ, t) for diffusion of particles in the energy space. For the simplified case of energy-independent transport coefficients, the nonlinear diffusion equation for the occupation-number distribution n ± (ǫ, t) becomes where the + sign represents bosons, and the − sign fermions. A stationary solution is given by the thermal distribution Eq. (2) for bosons and, correspondingly, the Fermi-Dirac distribution for fermions (quarks). In spite of its simple structure, the nonlinear diffusion equation with constant transport coefficients preserves the essential features of Bose-Einstein/Fermi-Dirac statistics which are contained in the quantum Boltzmann equation.
The transport equation can be solved exactly for a given initial condition n ± i (ǫ) using the nonlinear transformation outlined in Ref. [4]. For gluons, the nonlinear boson diffusion equation (NBDE) is more difficult to solve analytically due to the singularity at the chemical potential ǫ = µ < 0, and the need to consider the boundary conditions at the singularity. For fermions, there is no corresponding singularity, such that the exact solution of the nonlinear problem can be obtained with the free Green's function as performed in Refs. [4,12]. For gluons, the bounded solution n + (ǫ, t) [26] is briefly reviewed. It can be written as with the time-dependent partition function Z + (ǫ, t) ≡ Z(ǫ, t) obeying a linear diffusion equation In the absence of boundary conditions, the free partition function becomes The energy-independent prefactor a(t) in the partition function drops out when taking the logarithmic derivative in the calculation of the occupation-number distribution. The function F (x) with the integral over the initial conditions covers the full energy region −∞ < x < ∞.
Due to the occurence of the singularity, however, one eventually will have to consider boundary conditions at ǫ = µ < 0. For a solution without boundary conditions as in Refs. [4], Green's function G free (ǫ, x, t) of Eq. (6) is a single Gaussian but it becomes more involved once boundary conditions are considered. The function F (x) depends on the initial occupation-number distribution n i according to As discussed in Ref. [27], the definite integral can be replaced w.l.o.g. by the indefinite integral A i (x) over the initial distribution with ∂ x A i (x) = n i (y), resulting in For any given initial distribution n i , one can now compute the partition function and the overall solution for the occupation-number distribution function Eq. (5) analytically. The solution technique has been developed in Refs. [27,28] for the case of a cold bosonic atom gas that undergoes evaporative cooling. Here, and in Ref. [26] for different initial conditions, the approach is carried over to equilibrating gluons at relativistic energies. To solve the problem exactly, the chemical potential is treated as a fixed parameter. With lim ǫ↓µ n(ǫ, t) = ∞ ∀ t, one obtains Z(µ, t) = 0, and the energy range is restricted to ǫ ≥ µ. This requires a new Green's function that equals zero at ǫ = µ ∀ t. It can be written as (11) and the partition function with this boundary condition becomes The function F remains unaltered with respect to Eq. (10), except for a shift of its argument by the chemical potential. With a given initial nonequilibrium distribution n i , the NBDE can now be solved including boundary conditions at the singularity. The solution is given by Eq. (5).

A. Thermalization of Gluons
For massless gluons at the onset of a relativistic hadronic collision, an initial-momentum distribution n i (|p|) ≡ n i (p) = n i (ǫ) has been proposed by Mueller [29] based on Ref. [30]. It accounts, in particular, for the situation at the start of a relativistic heavy-ion collision [9]. It amounts to assuming that all gluons up to a limiting momentum Q s are freed on a short time scale τ 0 ∼ Q −1 s , whereas all gluons beyond Q s are not freed. Thus the initial gluon-mode occupation in a volume V is taken to be a constant up to Q s , as in Eq. (1) that was already used before in the relaxation-time approximation. Typical gluon saturation momenta for a longitudinal momentum fraction carried by the gluon x ≃ 0.01 turn out to be of the order Q s ≃ 1 GeV [13], which is chosen for the present model investigation.
Results for the gluon thermalization from n i (ǫ) to n eq (ǫ) according to Eq. (4) have been calculated in Ref. [4] for the free case, without considering boundary conditions at the singularity. As a consequence, diffusion into the negative-energy region occured, depleting the occupation in the infrared such that the asymptotic distribution differed from Bose-Einstein.
As a remedy, one has to extend the energy scale in Eq. (1) to µ ≤ ǫ < ∞, and include the boundary conditions at the singularity ǫ = µ < 0. This will cause the time-dependent solutions of the NBDE to properly approach the thermal Bose-Einstein distribution over the full energy scale as t → ∞. The initial condition is thus modified to include a δ-function singularity at ǫ = µ < 0 according to The δ-function singularity in the initial conditions of the NBDE has an analogous role as the singularity that can be added to the Boltzmann equation in order to act as a seed condensate [31] since the time evolution of the solutions without singularity does not lead to condensate formation. In Ref. [26], another choice of the initial conditions had been explored, with a thermal distribution in the negative-energy region that also has a singularity at ǫ = µ. Although the results differ in detail, the overall representation of thermalization is similar to the present results.
The time-dependent partition function with the above initial condition can now be calculated using the bound Green's function Eq. (11), and the function F (x) from Eq. (10). The latter contains an indefinite integral over the initial condition Eq. (13) that can be carried out to obtain (with x → x + µ in the argument of F (x) as required by the boundary conditions) The function F (x) is plotted in Fig. 2. Due to the singularity in its argument, F (x) has a discontinuity at x = 0. It is continuous, but not differentiable at x = Q s − µ. Both properties are essential to account for the equilibration near the singularity, and in the UV region. The Green's function of Eq. (11) that includes the IR boundary condition can explicitly be written as (15) With F (x) and G(ǫ, x, t), the partition function Z(ǫ, t) of Eq. (12) and its derivative ∂Z/∂ǫ can now be calculated, as well as the occupation-number distribution n(ǫ, t) from Eq. (5). The full calculation may in principle be carried out analytically. In the case of initial conditions that are appropriate for evaporative cooling of atomic Bose gases at very low energy, we have performed such an exact calculation including the boundary conditions at the singularity in Ref. [27]. Here I compute the partition function and its derivative using the NIntegrate and Derivative routines of Mathematica.

B. Discussion of the Solutions for Gluons
The bosonic equilibration time τ eq is taken as τ eq = 4D/(9v 2 ) ≃ 0.1 fm/c. This expression has been determined in Ref. [4] for a θ-function initial distribution in the UV. The approach to equilibrium provided by the solutions of the NBDE for the gluon distribution functions is shown in Fig. 3 at t = 6 × 10 −5 , 6 × 10 −4 , 6 × 10 −3 , 4 × 10 −2 , 0.12, and 0.36 fm/c, with decreasing dash lenghts. The steep cutoff in the UV at ǫ = Q s is smeared out at short times -this was the case already in the free solution without boundary conditions [4]. The diffusion coefficient is D = 1.17 GeV 2 c/fm, the drift coefficient v = −2.28 GeVc/fm. Correspondingly, x → x + µ as required by the boundary conditions, and a singularity at the origin. F (x) contains the integral over an initial nonequilibrium gluon distribution ni(x) according to Eq. (14). The parameters are given in the text.
the equilibrium temperature in this model calculation is T = −D/v ≃ 513 MeV = 0.513 GeV, as expected for the initial central temperature in a Pb-Pb collision at the LHC energy of √ s N N = 5 TeV [3]. Solutions for related, but different initial conditions (a thermal distribution in the energy region µ ≤ ǫ ≤ 0) have been discussed in Ref. [26].
The assumption of a constant negative chemical potential µ < 0 used in this work is, of course, an idealization that facilitates analytical solutions of the nonlinear problem. Here, the value of µ is calculated from particlenumber conservation with the density of states g(ǫ). For constant density of states and the above parameter values (T = 513 MeV), the result is µ = −0.08 GeV. For the more realistic density of states g(ǫ) ∝ ǫ 2 that is valid for a zero-mass relativistic system of gluons, the particle number is then also approximately conserved.
In general, particle-number conservation is strictly fulfilled e. g. for atomic Bose gases at much lower energy, but not for gluons in high-energy collisions. Driven by particle-number conservation, cold bosonic atoms can move into the condensed phase, thus diminishing the number of particles in the thermal cloud. The chemical potential in the equilibrium solution of the NBDE then becomes time dependent, as has been discussed in Ref. [27], albeit without a full quantum treatment of the condensed phase. It would become zero only in the limit of an inifinite number of particles in the condensed phase. Instead, it approaches a small but finite negative value for a finite number of particles.
In case of relativistic heavy-ion collisions, however, gluons can be created and destroyed. It is therefore unlikely that a condensed phase is actually formed, as had been proposed in model investigations where only soft elastic, number-conserving gluon collisions were considered [9]. Gluon condensate formation in relativistic collisions is essentially prevented by number-changing inelastic processes that correspond to splitting and merging of gluons, although a transient condensate formation is still being debated [10].
Hence, since inelastic collisions cannot be neglected, the gluon equilibrium distribution is expected to have a nearly vanishing, but still slightly negative, chemical potential which should be approached by the timedependent solutions of the NBDE. It would therefore be of interest to repeat the present calculation for a timedependent chemical potential, with µ(t) → 0 for t → ∞, as was done in Ref. [27] for the case of cold atoms. This requires, however, numerical work that goes substantially beyond the present analytic approach.

C. Discussion of the Solutions for Quarks
The local thermalization of valence quarks towards the Fermi-Dirac equilibrium distribution is easier to calculate in the nonlinear model because no singularity occurs for fermions. Hence, the free solutions can be used as was already discussed in Ref. [4] and in more detail in Ref. [12]. Hadron production in the early thermalization phase was implicitly considered, because the negative-energy region corresponds to particleantiparticle production. A typical result for the fermionic time-dependent occupation-number distribution as function of the transverse energy taken from Ref. [32] is shown in Fig. 4. Local thermalization for quarks occurs more slowly as for gluons due to Pauli's principle, and also because of the larger color factor for gluons.
The bounded solutions of the NBDE and the corresponding nonlinear fermion diffusion equation are, in par-  Figure 4. Occupation-number distribution n − (ǫ T, t) of a relativistic fermion system (valence quarks) as function of the transverse energy ǫ T including antiparticle production (ǫ T < 0). It is evaluated analytically from Eq. (4) at different times for the initial distribution ticular, tailored to local thermalization processes that occur in relativistic heavy-ion collisions at energies reached at RHIC and LHC. In the present example, they are applied to the local equilibration of quarks and gluons in central Pb-Pb collisions at a center-of-mass energy of 5 TeV per nucleon pair, leading to rapid thermalization with a local temperature of T ≃ 500 MeV. Since the thermalization occurs very fast -before anisotropic expansion fully sets in -, the analytic solution of the Figure 5. Schematic representation of the three-source model for relativistic heavy-ion collisions at RHIC and LHC energies in the center-of-mass system: Following the collision and slowing down (stopping) of the two Lorentz-contracted slabs (blue), the fireball region (center, yellow) expands anisotropically in longitudinal and transverse direction. At midrapidity, it represents the main source of particle production. The two fragmentation sources (red) contribute to particle production, albeit mostly in the forward and backward rapidity regions.
In stopping, these are the only sources. From Ref. [33].
problem in 1+1 dimensions appears permissible. The hot system will subsequently expand anisotropically and cool rapidly, as is often modeled successfully by relativistic hydrodynamics [1], until hadronization is reached at T ≃ 160 MeV. Further refinements of the thermalization model such as time-dependent transport coefficients are conceivable, but are unlikely to allow for analytic solutions. A microscopic calculation of the transport coefficients with an investigation of their dependencies on energy and time would be very valuable. Extensions of the NBDE itself to higher dimensions in order to account for possible anisotropies should also be investigated.

III. STOPPING: NET-PROTON DISTRIBUTIONS
The incoming baryons with energies available at SPS, RHIC or LHC are being stopped on an equally short time scale as the local equilibration occurs: In the course of the collision shown schematically in Fig. 5, the system is being slowed down, essentially through collisions of the incoming valence quarks with soft gluons in the respective other nucleus. Various models to account for this process and its energy dependence have been developed, in particular in Refs. [13,14] and related works which are relying on the appropriate parton distribution functions (PDFs). They yield agreement with the available netproton (proton minus antiproton) stopping data. Different from the nonequilibrium-statistical approach to initial thermalization, such models do not consider a time with NA49 data [34]. Red curves are the initial distributions broadened by the Fermi momentum, the final distribution is from a QCD-inspired model [13] with the saturation-scale exponent λ = 0.2 and Q 2 0 = 0.09 GeV 2 (see text). It agrees with the data, and corresponds to the distribution at the interaction time t = τint in a time-dependent formulation. The six intermediate solid curves at t/τint = 0.01, 0.02, 0.05, 0.1, 0.2 and 0.5 account for the time dependence in a nonequilibriumstatistical relativistic theory. From Hoelck and Wolschin [35].
dependence. However, by using the rapidity distribution calculated from the PDFs, and the initial valence-quark distribution, one can, in addition, account for the time development from the initial to the final distribution with an appropriate fluctuation-dissipation relation.
The fragmentation peaks in stopping occur mainly due to the interaction of valence quarks with soft gluons in the respective other nucleus. For net protons (protons minus antiprotons), their positions in rapidity space y = 0.5 ln[(E || +p || )/(E || −p || )] with m p = m p can be obtained from Ref. [13], and references therein, as This expression accounts for the peak in the forward region, and there is a corresponding one with y → −y for the backward peak. The longitudinal momentum fraction of the valence quark v that experiences stopping is , the one for the soft gluon g in the target is . The distribution function of the valence quarks is q v (x 1 , p T ), the one of the gluons is f g (x 2 , p T ). The latter represents the Fourier transform of the forward dipole scattering amplitude N (x 2 , r T ) for a quark dipole of transverse size r T [14]. To account for the correct normalization in net-proton or net-baryon distributions, the normalization constant C is adjusted such that the integral of Eq. (17) agrees with the total number of participant protons or baryons.
Fragmentation peaks are then found to occur at y = ±y peak in rapidity space. At suffiently high energyin particular, at LHC energies -these positions become sensitive to the gluon saturation scale Here, A is the mass number, Q 0 the momentum scale, x < 1 the momentum fraction carried by the gluon, and λ the saturation-scale exponent. values of the gluon saturation scale and found to agree with net-proton data from SPS and RHIC. We are presently incorporating these results into a timedependent nonequilibrium-statistical relativistic theory [35]. A typical outcome of our approach for central Pb-Pb at SPS energies of √ s N N = 17.3 GeV is shown in Figure 6 where it is compared with NA49 data [34] in rapidity space. The initial distributions are shown as red curves, with the initial broadening due to the Fermi motion. The final distribution agrees with the one from our QCD-inspired model [13] with a saturation-scale exponent λ = 0.2 and Q 2 0 = 0.09 GeV 2 , translating into a gluon saturation momentum of Q s ≃ 0.55 GeV/c at x = 10 −4 . At the interaction time t = τ int , the timedependent distribution agrees with the data. The six intermediate solid curves at t/τ int = 0.01, 0.02, 0.05, 0.1, 0.2 and 0.5 show the time dependence in our nonequilibriumstatistical relativistic theory [35].
A larger gluon saturation momentum Q s was found to produce more stopping, as does a larger mass number A [13]. In the context of an investigation of particle production, the agreement between the calculated stopping distributions and the data will be taken as evidence for the importance of fragmentation contributions also in charged-hadron production.
In Ref. [36], Mehtar-Tani and I found that the fragmentation peak positions in stopping depend in a large c.m. energy range 6.3 GeV ≤ √ s N N ≤ 200 GeV linearly on the beam rapidity y beam and the saturation-scale exponent λ according to At the current LHC energy of 5.02 TeV Pb-Pb corresponding to y beam = ± ln( √ s N N /m p ) = ± 8.586 and with a gluon saturation-scale exponent λ ∼ 0.2 one therefore expects y peak ≃ ± 6. Due to the lack of a suitable forward spectrometer at LHC, the rapidity region of the peaks will thus not be accessible for identified protons in the coming years at LHC energies. Nevertheless, the partonic processes that mediate stopping also contribute to hadron production at LHC energies and hence, one expects fragmentation events in particle production. Experience with stopping versus hadron production at SPS and RHIC energies [15,16] has shown that the fragmentation peaks in particle production occur consistently at somewhat smaller absolute rapidities than the ones in stopping.
In net-baryon (proton) distributions charged baryons produced from the gluonic source cancel out because particles and antiparticles are generated in equal amounts. In charged-hadron production, however, this is not the case. Instead, three sources contribute provided the energy is sufficiently high, √ s N N > 20 GeV. It was found in Ref. [16] that the dependence of their parti-cle content on c.m. energy differs: The fragmentation sources contain N qg ch ∝ ln(s N N /s 0 ) charged hadrons, whereas the midrapidity-centered source that arises essentially from the interaction of low-x gluons contains N gg ch ∝ ln 3 (s N N /s 0 ) charged hadrons. Due to the strong ∝ ln 3 dependence, it becomes more important than the fragmentation sources at LHC energies [37].
With the three sources, the total rapidity distribution for produced charged hadrons becomes The fragmentation distributions R 1,2 (y, t) and gluonic distributions R gg (y, t) can be calculated in a timedependent phenomenological model such as the relativistic diffusion model [15], or in microscopic theories. The strong interaction ceases to act at the interaction time (≡ freezeout-time) t = τ int and theoretical distributions may be compared to data in a χ 2 -optimization. A transparent phenomenological model to calculate and predict the distribution functions of produced particles is the relativistic diffusion model [15,38]. In the RDM, the initial distribution functions are evolved up to τ int /τ y with the rapidity relaxation time τ y using the analytical moments equations. The mean values y 1,2 of the fragmentation distributions that are related to τ int /τ y can be determined from the data. The details of the model calculations are given in the corresponding Refs. [15,16,37].

IV. CHARGED-HADRON PRODUCTION
A. Transverse-momentum distributions As discussed frequently in experimental papers [39] and also in Ref. [16], the transverse-momentum distributions of produced charged hadrons in relativistic heavyion collisions show an exponential behaviour in the thermal regime as accounted for by the Maxwell-Jüttner distribution [40] with the modified Bessel function of the second kind K 2 (m/T ), the Lorentz-factor freeze-out temperature T and hadron mass m. Beyond p T ≃ 4 GeV/c, however, a transition to a power-law (straight lines in a log-log plot) occurs. It is attributed mostly to the recombination of soft partons, and fragmentation of hard partons. In addition to detailed theoretical approaches, this transition can be modelled phe-nomenologically using distribution functions of the QCDinspired form (see references in Wilk and Wong [41]) used by Hagedorn [42] for high-energy pp and pp collisions with a normalization constant C and parameters p 0 , n. This expression describes the transition from exponential for p T → 0 as in the Jüttner distribution Eq. (21) (with p 0 = nT and p T → m T ), to power-law behaviour (∝ (p T /p 0 ) −n for p T → ∞) . Using Eq. (23), Fig. 7 shows p T -distributions of produced charged hadrons at four centralities in 5.02 TeV Pb-Pb compared with ALICE data [43] (peripheral spectra are scaled for better visibility; statistical and systematic error bars are smaller than the symbol size). The data are well represented through many orders of magnitude with a power index n = 8.2 and p 0 = 3 GeV/c (Fig. 7), but at high p T deviations occur which are attributed to hard processes that require a pQCD treatment. This corresponds to the occurence of a minimum in the nuclear modification factor for produced charged hadrons as function of p T found already at √ s N N = 2.76 TeV [44], and confirmed at 5.02 TeV [43].
There is presently no theoretical derivation for the value of the power index n that is needed to reproduce the experimental p T -distributions in relativistic heavyion collisions. It is therefore not obvious from the present analysis which fraction of low-p T particles could be due to nonequilibrium processes that differ from thermal emission out of a single expanding fireball. In particular, one can not distinguish particles emitted from the fireball and those arising from the fragmentation sources at low p T . Hence, the analysis of transverse momentum distributions in terms of Eq. (23) is presently only suitable to distinguish high-p T hard events from the bulk of (thermal and nonequilibrium) charged-hadron emission.

B. Pseudorapidity distributions
The distinction of particles emitted from the fireball and those from the fragmentation sources is more transparent in rapidity or pseudorapidity distributions of produced charged hadrons. The existence of the fragmentation sources is evident from the measurements of stopping in heavy-ion collisions as discussed in Section 3. Pseudorapidity distributions dN ch /dη with η = − ln [tan(θ/2)] depend only on the scattering angle θ. Particle identification is not needed here and hence, they are much easier to obtain at large η-values (small scattering angles) compared to rapidity distributions.
The pseudorapidity distributions dN ch /dη for produced charged hadrons in relativistic heavy-ion collisions emerge from a superposition of the fragmentation sources and a midrapidity source. According to the discussion in Ref. [37] and at the end of Section 3, the particle content of the low-x gluon source rises rapidly according to N gg ∝ ln 3 (s N N /s 0 ). To convert rapidity distributions dN ch /dy to pseudorapidity distributions dN ch /dη, the corresponding Jacobian is required The hadron mass is m and the transverse momentum p T . Since the transformation depends on the squared ratio (m/p T ) 2 of mass and transverse momentum of the produced particles, its effect increases with the mass of the particles and is most pronounced at small transverse momenta. In Ref. [45], we have determined the Jacobian J 0 at η = y = 0 in central 2.76 TeV Pb-Pb collisions for identified π − , K − , and antiprotons from the experimental values dN dη | exp and dN dy | exp as J 0 = 0.856. We can solve Eq. (25) for p T ≡ p eff T to obtain The mean mass m can be calculated from the abundancies of pions, kaons, and antiprotons. Using J 0 , the Jacobian can be written independently from the values of m and p eff T as J (η, J 0 ) = cosh(η) The result for central 2.76 TeV Pb-Pb collisions was found in Ref. [45] to be J(η) = cosh(η)[1.365 + sinh 2 (η)] −1/2 .
The Jacobian affects the pseudorapidity distributions mostly near midrapidity, where it generates a dip as is obvious from Figure 8: A prediction in the RDM with linear drift from Ref. [16] is compared with ALICE data for central Pb-Pb at 5.02 TeV [46] (dashed upper curve), and a χ 2 -optimization within the five parameter RDM is carried out, upper solid curve. It differs only slightly from the prediction.
The nonequilibrium evolution of all three partial distribution functions R k (y, t) (k = 1, 2, gg) towards the thermodynamic equilibrium (Maxwell-Jüttner) distribu- is accounted for in the Relativistic Diffusion Model [15,16,38,47] through solutions of the Fokker-Planck equation The drift functions J k (y, t) and diffusion functions D k (y, t) depend on rapidity and time. However, if the diffusion coefficients are taken as constants D k , and the drift functions assumed to be linearly dependent on the rapidity variable y, the FPE acquires the Ornstein-Uhlenbeck form [48] which can be solved analytically in rapidity space [38]. In this case it is easy to show that for t → ∞ all three subdistributions approach a single Gaussian in y space which is centered at midrapidity y = 0 for symmetric systems, or at the appropriate equilibrium value y = y eq for asymmetric systems [15,16,37]. For stopping, only the two fragmentation distributions contribute, approaching the thermal equilibrium distribution for t → ∞, as discussed in the previous section. If the system reaches a stationary distribution that differs from the thermal one as discussed in Section 3 for the case of stopping calculated in a QCD-based model, the underlying fluctuation-dissipation relation becomes more complicated [35]. To derive the above FPE (29) in the context of relativistic heavy-ion collisions, one can make use of a theory for a special class of non-markovian processes in spacetime discussed in Refs. [49,50], which are equivalent to relativistic Markov processes in phase space (RMPP). These processes give rise to a generalised FPE which is suitable for describing relativistic diffusive particle dynamics. Our Eq. (29) is conceptually a special case of such a RMPP formalism for charged-hadron production in rapidity space. An application of RMPP to baryon stopping will be shown in Ref. [35].
The t → ∞ limit of the FPE solution for constant diffusion and linear drift is found to deviate slightly from the Maxwell-Jüttner distribution. The discrepancies are small and become visible only for sufficiently large times. To ensure that the asymptotic solution yields the Maxwell-Jüttner distribution Eq. (28), a RDM with the sinh-drift is required as was discussed in Refs. [51,52]. The corresponding fluctuation-dissipation relation (FDR) that connects drift and diffusion becomes [52] A k = m T D k /T .
If the asymptotic distribution is not Maxwell-Jüttner, but -as in the case of stopping, where it may be provided by a QCD-based distribution from a calculation as performed in Ref. [13] -, a different form of the fluctuationdissipation relation will result, as will be discussed in Ref. [35]. The strength of the drift force in the fragmentation sources k = 1, 2 depends on the distance in y-space from the beam rapidity, which enters through the initial conditions. With Eq. (28), the rapidity distribution at thermal equilibrium can then be derived [53] as where C is proportional to the overall number of produced charged hadrons N tot ch , or -in case of stoppingto the number of net baryons (protons) in the where C is proportional to the overall number of produced charged hadrons N tot ch , or -in case of stopping -to the number of net baryons (protons) in the respective centrality bin. Since the actual distribution functions remain far from thermal equilibrium, the total particle number is evaluated based on the nonequilibrium solutions of the FPE, which are adjusted to the data in χ 2 -optimizations.
In particular, one can determine the drift amplitudes A k from the position of the fragmentation peaks as inferred from the data, and then calculate theoretical diffusion coefficients as D k = A k T /m T . Since the fireball and both fragmentation sources also expand collectively, the actual distribution functions will, however, be broader than what is obtained from Eq. (31). To account for this broadening through collective expansion, we use diffusion coefficients (or widths of the partial distributions) that are adapted to the data in both stopping and particle production. From the integral of the overall distribution function the total particle number can be obtained. In case of the FPE with sinh-drift, the FPE must be solved numerically as described in Ref. [52].
Results for central collisions of symmetric systems in the three-source RDM with linear drift are summarized in Figure 8. The dependence of the pseudorapidity distributions on c.m. energy in central Au-Au collisions at 19.6 GeV, 130 GeV, and 200 GeV RHIC energies as well as in central Pb-Pb at 2.76 TeV and 5.02 TeV LHC energies are displayed. In addition to RDM calculations with parameters for the lower energies from [15] compared with data from Refs. [18,44,54], a prediction for 5.02 TeV Pb-Pb (dashed upper curve) from Ref. [16] is compared with recent ALICE data at 0-5 % centrality [46]. A fit to the data within the five-parameter RDM is also displayed (solid curve).
Only the fragmentation sources contribute at the lowest RHIC energy of 19.6 GeV that is shown here -which is comparable to the highest SPS energy of 17.3 GeV -(see dot-dashed curves), but at higher energies the gluonic source rapidly rises and becomes the largest source of particle production at an energy of ∼ 2 TeV, which is between energies reached at RHIC and LHC.
I have investigated the dependence of the particle content of the three sources on center-of-mass energy per particle pair √ s N N in Ref. [37]. The gluonic source is absent for √ s N N 20 GeV, see the 19.6 GeV Au-Au PHO-BOS result in Figure 8. At this relatively low energy, charged-hadron production arises only from the fragmentation sources which overlap in rapidity space and hence, appear like a single gaussian ("thermal") source. The total charged-hadron production has been found experimentally to depend linearly on ln(s N N /s 0 ), see for example central Pb-Pb NA50 data at 8.7 GeV and 17.3 GeV [55] and low-energy Au-Au PHOBOS results [18].
In my RDM-analysis with three sources published in Ref. [37], it has turned out that the dependence of the fragmentation sources N qg ch ∝ ln(s N N /s 0 ) indeed continues at higher energies, see Figure 9. In addition to what was shown in Figure 4 of Ref. [37] with an extrapolation to 5.02 TeV, now my results based on an analysis of new ALICE 5.02 TeV data [57] are included in this figure. The gluonic source is confirmed to have a strong energy dependence N gg ch ∝ ln 3 (s N N /s 0 ). As discussed in Ref. [37], the rise of the cross section in the central distribution is driven by the growth of the gluon density at small x and theoretical arguments [59] suggest a ln 2 s asymptotic behaviour that satisfies the Froissart bound [60]. Because the beam rapidity is ∝ ln(s N N ), the integrated yield from the gluonic source then becomes proportional to ln 3 s, in agreement with the phenomenological analysis, and the new 5.02 TeV data. It was mentioned in Ref. [37] that there exist also further experimental confirmations of this result at RHIC energies based on STAR data for dijet production, see [61] and references therein.    Figure 9. The total charged-hadron production in central Au-Au and Pb-Pb collision in the energy region 19.6 GeV to 5.02 TeV is following a power law (solid upper line), whereas the particle content in the fragmentation sources is Nqg ∝ ln (sNN /s0), dash-dotted curve. The particle content in the mid-rapidity source obeys Ngg ∝ ln 3 (sNN /s0), dashed curve. The energy dependence of the measured mid-rapidity yields is shown as a dotted line, with PHOBOS data [54] at RHIC energies, and ALICE data [57,58]  The 5.02 TeV Pb-Pb data confirm that the sum of produced charged hadrons integrated over η is close to a power law N tot ch ∝ (s N N /s 0 ) 0.23 with s 0 = 1 TeV 2 as shown in Figure 9 for central Au-Au and Pb-Pb col-lisions, upper line. At RHIC energies Busza noticed that the integrated charged-particle multiplicities scale as ln 2 (s N N /s 0 ) [62], but the energy dependence up to the maximum LHC energy of 5.02 TeV has turned to be even stronger due to the high gluon density. In Ref. [37] it was shown that the midrapidity yields for central Au-Au and Pb-Pb collisions are dN tot ch dη η≃0 = 1.15 · 10 3 (s N N /s 0 ) 0.165 (33) with s 0 = 1 TeV 2 (dotted line, data points from Phobos [54] and ALICE [57,58]). More detailed aspects of the interplay between fragmentation sources and gluonic source appear when investigating the centrality dependence of chargedhadron pseudorapidity distributions, as has been done in Refs. [63,64] for the asymmetric systems 200 GeV d-Au and 5.02 TeV p-Pb, and in Refs. [15,45] for 2.76 TeV Pb-Pb.

C. Limiting Fragmentation at RHIC and LHC energies
Using the RDM with both, linear and sinh-drift, we have investigated whether the limiting fragmentation conjecture is fulfilled at energies reached at RHIC and LHC [33]. The significance of the fragmentation region in relativistic heavy-ion collisions had been realized when data on Au-Au collisions in the energy range √ s N N =19.6 GeV to 200 GeV became available at RHIC [17][18][19]. The pseudorapidity distributions of produced charged particles for a given centrality bin scale with energy according to the limiting fragmentation (LF), or extended longitudinal scaling, hypothesis: Over a large range of pseudorapidities η ′ = η − y beam in the fragmentation region with the beam rapidity y beam , the chargedparticle pseudorapidity distribution is found to be energy independent.
The phenomenon was first shown to be present in pp data, in a range from 53 up to 900 GeV [65], following a prediction for hadron-hadron and electron-proton collisions in Ref. [66]. With increasing collision energy the fragmentation region grows in pseudorapidity space. It can cover more than half of the pseudorapidity range over which particle production occurs. Especially in relativistic heavy-ion collisions, the approach to a universal limiting curve is a remarkable feature of the particle production process.
As discussed in Ref. [33] and references cited therein, it is an interesting question whether limiting fragmentation will persist at the much higher incident energies that are available at the LHC, namely, √ s NN = 2.76 and 5.02 TeV in Pb-Pb collisions. At these energies, experimental results in the fragmentation region are not available due to the lack of a dedicated forward spectrometer. If one wants to account for the collision dynamics  Figure 10.
Comparison of the three-source RDMdistributions with linear and sinh-drift, PHOBOS data [54], and ALICE data [46,56]  The zoom into this region shows that the RDM with sinh-drift is consistent with limiting fragmentation at RHIC and LHC energies. From Kellers and Wolschin [33], where details and parameters are given.
more completely, however, this region is most interesting. We have therefore investigated in Ref. [33] to what extent limiting fragmentation can be expected to occur in heavy-ion collisions at LHC energies. The result from that investigation is summarized in Figure 10: Limitingfragmentation scaling can be expected to hold at both, RHIC and LHC energies. This conclusion agrees with microscopic numerical models such as AMPT [67], but it disagrees with expectations from simple parametrizations of the rapidity distributions such as the difference of two Gaussians. It also disagrees with predictions from the thermal model, which does not explicitly treat the fragmentation sources but refers only to particles produced from the hot fireball. In contrast, in our approach the fragmentation sources play an essential role. Only future upgrades of the detectors would make it possible to actually test the limiting-fragmentation conjecture experimentally at LHC energies.

V. QUARKONIA AND THE QGP
Among the hard probes in relativistic heavy-ion collisions, the modification of quarkonia yields in the presence of the quark-gluon plasma has an outstanding role. Charmonia (J/ψ) suppression due to the screening of the real part of the Cornell-type quark-antiquark potential in the hot medium had initially been suggested by Matsui and Satz in 1986 as a QGP signature [23]. It was realized later that the potential in the medium is an optical one,  Figure 11. Transverse-momentum dependence of the suppression factor RAA(Υ(1S)) for the spin-triplet ground state in minimum-bias Pb-Pb collisions at √ sNN = 5.02 TeV.
The (upper) dashed curve shows the suppression in the hot medium, the (lower) solid curve the suppression including reduced feed-down, which is important for the ground state, but not for excited states. The theoretical prediction is from Ref. [3], data are from CMS [69].
with the imaginary part [24] causing dissociation of the quarkonia states in the hot medium and thus, quarkonia suppression when compared to the production rates from pp collisions at the same center-of-mass energy scaled with the number of binary collisions. In addition to the associated collisional damping widths of the quarkonia states, thermal gluons can dissociate these states, and the gluon-induced dissociation widths [25] can be treated separately from the damping [68]. An important role is played by the reduction of feed-down in the heavyion case as compared to pp, because feed-down from the higher to the lower states is hindered if the higher states are screened away, or depopulated. If the medium contains a large number of heavy quarks as is the case for charm quarks in Pb-Pb collisions at LHC energies, statistical recombination cannot be neglected at sufficiently low transverse momentum, and there is an interplay of dissociation and recombination. There is meanwhile a large number of publications and reviews about charmonia physics in relativistic heavy-ion collisions available [39]. Bottom quarks are about three times heavier than charm quarks, have a correspondingly smaller production cross section, and hence, are less abundant even at LHC energies. Statistical recombination is therefore less important, and expansions in terms of (1/m) are more precise. Consequently, bottomonia provide a cleaner probe of the QGP properties such as the initial central temperature, and have been investigated in detail both theoretically and experimentally. The bottom quarks are produced on a very short time scale of 0.02 fm/c in the initial stages of the collision, before the QGP of light quarks and gluons is actually being produced. The formation time of bottomonia states is larger, in the range τ F ≃ 0.3 − 0.6 fm/c. It is less precisely known, may differ for the individual states such as Υ(1S, 2S, 3S) and Theoretical predictions are from Hoelck, Nendzig, Wolschin [3], CMS data are from Ref. [69].
χ b (1P, 2P, 3P ), and could depend on the temperature of the emerging QGP, which would further enlarge it [70].
Since the spin-triplet Υ(1S) state is particularly stable with a binding energy of ≃ 1.1 GeV, it has a sizeable probability to survive as a color-neutral state in the colored hot quark-gluon medium of light quarks and gluons that is created in a central Pb-Pb collision at LHC energies, even at initial medium temperatures of the order of 400 MeV or above.
There exists a considerable literature on the dissociation of quarkonia, in particular of the Υ meson [71][72][73], in the hot quark-gluon medium; see Ref. [74] and references therein for a review. In minimum-bias Pb-Pb-collisions at LHC energies of √ s N N = 5.02 TeV in the midrapidity range, the strongly bound Υ(1S)-state is found to be suppressed down to about 38 % as compared to the expectation from scaled pp collisions at the same energy. The Υ(1P ) state has a smaller binding energy and is even more suppressed, down to 12 % [69].
Various theoretical approaches such as Refs. [75][76][77][78], and their recent updates to higher energy, are available that allow for an interpretation of the data. In the next section, results of our model are reported that aims to account for the previous Pb-Pb results at 2.76 TeV and to predict results for the higher energy of 5.02 TeV, which are then compared to the data that have meanwhile become available.
A. Υ(1S, 2S) suppression in Pb-Pb at LHC energies In Refs. [3,68,79] we have devised a model that accounts for the screening of the real part of the potential, the gluon-induced dissociation of the various bottomonium states in the hot medium (gluodissociation), and the damping of the quark-antiquark binding due to the presence of the medium which generates an imaginary part of the temperature-dependent potential. Screening is less important for the strongly bound Υ(1S) ground state, but it is relevant for the bb excited states, and also for all cc bound states.
Due to screening and depopulation of the excited states in the hot medium, the subsequent feed-down cascade towards the Υ(1S) ground state differs considerably from what is known based on pp collisions. The LHCb collaboration has measured a feed-down fraction of Υ(1S) originating from χ b (1P ) decays in pp collisions at √ s = 7 TeV of 20.7 % [80], and the total feed-down from excited states to the ground state is estimated to be around 40 % [81] at LHC energies. If feed-down was completely absent because of screening and depopulation of excited states in the hot medium, a modification factor R AA (Υ(1S)) ≃ 0.6 would thus result, whereas the measured modification factor of the Υ(1S) state in minimum-bias Pb-Pb collisions at 2.76 TeV is R AA (Υ(1S)) = 0.453 ±0.014 (stat)±0.046 (syst) [82], and at 5.02 TeV R AA (Υ(1S)) = 0.378 ± 0.013 (stat)±0.035 (syst) [69]. Hence, there clearly exist in-medium suppression mechanisms for the strongly bound Υ(1S) state which we aim to account for in detail, together with the suppression of the excited states, and the reduced feed-down.
In our model calculation [3], we thus determine the respective contributions from in-medium suppression, and from reduced feed-down for the Υ(1S) ground state, and the Υ(2S) first excited state in Pb-Pb collisions at both LHC energies, 2.76 TeV and 5.02 TeV. The p Tdependence and the role of the relativistic Doppler effect on the measured transverse-momentum spectra are considered. For the Υ(2S) state, the QGP effects are expected to be much more important than reduced feeddown. We compare in Ref. [3] with centrality-dependent CMS data [71,82] for the Υ(1S) and Υ(2S) states in 2.76 TeV Pb-Pb collisions, and predict the p T -and centrality-dependent suppression at the higher LHC energy of √ s N N = 5.02 TeV. The predictions at 5.02 TeV are compared with recent CMS data in Fig. 11 for the transverse-momentum dependence, Fig. 12 for the centrality dependence of Υ(1S), and Fig. 13 for the centrality dependence of Υ(2S).
For symmetric systems such as Au-Au at RHIC or Pb-Pb at LHC, we do not include an explicit treatment of CNM effects such as shadowing in the present study. These are, however, important in asymmetric collisions such as p-Pb where most of the system remains cold during the interaction time, and we have considered them in our corresponding calculations [83]. Statistical recombination of the heavy quarks following bottomonia dissociation is disregarded: Although this is certainly a relevant process in the J/ψ case, the cross section for Υ production is significantly smaller.
The anisotropic expansion of the hot fireball is accounted for using hydrodynamics for a perfect fluid that includes transverse expansion. Such a simplified nonviscous treatment [3,68] of the bulk evolution appears to be tolerable because conclusions on the relative importance of the in-medium suppression versus reduced feeddown are not expected to depend much on the details of the background model. When calculating the in-medium dissociation, we consider the relativistic Doppler effect that arises due to the relative velocity of the bottomonia with respect to the expanding medium. It leads to more suppression at high p T , and to an overall rather flat dependence of R AA on p T .
Our predictions for the p T -dependent Υ-suppression in 5.02 TeV Pb-Pb collisions are shown together with recent CMS data [69] in Fig. 11; see the caption for details. The in-medium modification factor (dashed) first rises at small transverse momentum, because escape from the hot zone becomes more likely with increasing p T , but then falls off when the increase in effective temperature becomes more pronounced. For the Υ(1S) state, a substantial fraction of the suppression, in particular at low p T , is due to reduced feed-down, solid curve. The corresponding centrality-dependent suppression (integrated over p T ) is shown in Fig. 12 to be in agreement with the data [69] for the Υ(1S) state. Related ALICE data at more forward rapidities 2.5 < y < 4 are roughly consistent within the error bars [84]. The suppression of the Υ(2S) state in Fig. 13 is mostly in-medium, with only a small contribution due to reduced feed-down. The prediction shows less suppression than the data in peripheral collisions. We have shown in Ref. [85] that the extra suppression of the loosely bound Υ(2S) state is most likely not due to the strong electromagnetic fields in more peripheral collisions. Hence, the origin of this effect is presently unknown.
It is of considerable interest to determine if the bottomonia distributions in more peripheral collisions become anisotropic, as has been found for produced particles in general. The quadrupole part of the momentum anisotropy is due to the almond-shaped spatial anisotropy of the overlap region, which translates to momentum space. It is more pronounced for lighter mesons such as pions, and can be quantified by the ellipticity v 2 of the momentum distribution in a Fourier decomposition of the experimentally determined, event-averaged particle distribution [1] with the azimuthal angle φ, the mean flow angle ψ n , and N the mean number of particles of interest per event (charged hadrons or identified particles of a specific species). The flow planes are, however, not experimentally known, and hence, the anisotropic flow coefficients are obtained using azimuthal angular correlations between the observed particles. The experimentally reported anisotropic flow coefficients from two-particle correlations can then be obtained as the root-mean-square values, v n {2} ≡ v n ≡ v 2 n , and the flow coefficients are being measured not in individual events, but in centrality classes.
Whereas for charged hadrons the flow coefficients have been measured very precisely [1] with v 2 -values up to 20 − 30% for pions, kaons, and antiprotons and maxima near p T ≃ 3 GeV/c, this is more difficult for quarkonia due to the much smaller production rates. For the charmonium ground state, the ellipticity coefficient at forward rapidity (2.5 < y < 4) in the centrality class 5 − 60% is v 2 ≃ 3 − 8%, depending on p T [86], implying that J/ψ shows elliptic flow, albeit on a smaller scale due to the larger mass.
It may, therefore, appear possible that even bottomonium exhibits flow, following the mass ordering of lighter particles resulting from a collective expansion of the medium [90]. Indeed, the large statistical error bars on the presently available bottomonium data for v 2 [86] do not yet exclude such a possibility, although it is quite doubtful whether this meson -which is about three times heavier than charmonium -flows with the expanding hot medium. Instead, the Υ(1S) is expected to essentially maintain its trajectory in the hot QGP, unless it is dissociated. Still, its momentum distribution in more peripheral collisions may exhibit a finite v 2 due to   [88] in the backward (Pb-going, left) and forward (p-going, right) region, for minimum-bias centrality. Results for CNM effects that include shadowing, energy loss, and reduced feed-down (dashed curves, blue) are shown together with calculations that incorporate also QGP effects (solid curves, red). The error bands result from the uncertainties of the parton distribution functions that enter the model. Calculations from Ref. [83].
the anisotropic escape from the fireball, because the path length from Υ formation to escape in the transverse plane depends on the azimuthal angle. This mechanism had already been suggested for J/ψ by Wang and Yuan at SPS and RHIC energies [91], and has been used in Ref. [92] for the bottomonium states in non-central 2.76 TeV Pb-Pb. There, the maximum of v 2 (Υ(1S)) including feed-down contributions is found to be below 1% in the 40 − 50% centrality class. We have performed a corresponding calculation within our model in the 5 − 60% centrality class and compare the p T dependence with the available AL-ICE data [86] , see Figure 14. The anisotropy is very small, compatible with zero. For more definite conclusions, one has to wait for a reduction of the experimental uncertainties in run 3.

TeV
Regarding bottomonia in asymmetric collisions, p-Pb at √ s N N = 8.16 TeV has been investigated experimentally by the LHCb [88] and ALICE [93] collaborations, and cold nuclear matter predictions had been published by a group of theorists [94]. Clearly, CNM effects are much more relevant than in symmetric systems, because the bulk of the hadronic matter remains cold during the interaction, see Fig. 15. The most relevant CNM effect is the modification of the parton distribution functions in the nuclear medium, which have been studied by many authors. A typical result for the PDF modifications with shadowing at small values of Bjorken-x, and antishadowing at intermediate x-values as obtained with EPPS16 [87] is shown in Fig. 16. Shadowing causes a reduction of the Υ(nS) yields in p-Pb as compared to scaled pp, whereas antishadowing results in an enhancement. Shadowing is somewhat more pronounced if one, in addition, considers coherent energy-loss mechanisms in the cold medium. Still, these are not sufficient to interpret the available data in terms of CNM effects, as becomes obvious from direct comparisons, in particular, for the Υ(2S) state.
There is, however, a spatially small hot zone (fireball) with an initial central temperature that is comparable to the one in a symmetric system (Fig. 17), and during its expansion and cooling, it contributes to bottomonia dissociation in regions where the temperature remains above the critical value. We have investigated the respective cold-matter and hot-medium effects on Υdissociation in 8.16 TeV p-Pb collisions in Ref. [83]. Representative results from this work are shown in Fig. 18 for the transverse-momentum dependence, and Fig. 19 for the rapidity dependence. The plots show CNM (blue, upper bands) and CNM plus QGP (red, lower bands) effects  Figure 19. Calculated rapidity-dependent nuclear modification factors R p-Pb for the Υ(1S) (top) and Υ(2S) state (bottom) in p-Pb collisions at √ sNN = 8.16 TeV with preliminary ALICE data [89] (triangles) and LHCb data [88] (circles). Results for CNM effects that include shadowing, energy loss, and reduced feed-down (dashed curves, blue) are shown together with calculations that incorporate also QGP effects (solid curves, red). The error bands result from the uncertainties of the parton distribution functions that enter the calculations. From Dinh, Hoelck, Wolschin [83].
on the Υ(1S) and Υ(2S) yields in 8.16 TeV p-Pb collisions at the LHC. The transverse-momentum dependence in the backward direction (top) shows enhancement due to antishadowing when only the CNM effects are considered, whereas the data for Υ(1S) are clearly suppressed at p T < 10 GeV/c and for Υ(2S) at all measured p T values. This discrepancy is cured through the consideration of the momentum-dependent dissociation in the QGP, as shown in our cold-matter plus hot-medium calculation (red). The forward/backward asymmetric shape of the nuclear modification factors as functions of rapidity (Fig. 19) arises from the different cold-matter effects in the forward and backward regions (in particular, shadowing/antishadowing of the parton distribution functions, but also energy loss in the relatively cold medium). The additional suppression due to the dissociation in the hot fireball is again shown in the lower (red) curves, which are in better agreement with the data for the Υ(1S) ground state not only in the backward, but also in the forward direction. The substantial role of the hot-medium effects is even more pronounced for the Υ(2S) first excited state, where the CNM calculation shows enhancement in the backward region, whereas the full calculation with in-medium dissociation displays a suppression down to about 70% -in agreement with the LHCb data [88] and the ALICE data point [93].
There have been attempts to explain the discrepancy between CNM calculations and data for the Υ(nS) suppression in p-Pb in terms of interactions with comoving hadrons, in particular, pions [95]. We have not included this process in our calculations -initially on the grounds that interactions of the bottomonia states with comovers were found to be unimportant at LHC energies in the work of Ko et al. [96] about Υ absorption in hadronic matter. Probably one eventually has to consider both, comover interactions plus suppression in the hot QGP zone in order to fully understand the Υ modification data in asymmetric systems.

VI. CONCLUSIONS
This article presents aspects of relativistic heavy-ion collisions with an emphasis on energies reached at the Relativistic Heavy Ion Collider RHIC and the Large Hadron Collider LHC. It does not attempt to be a review of the field, which is available in recent textbooks and in the proceedings of Quark Matter conferences such as Ref. [39]. Instead, a specific phenomenological viewpoint with experimental data as a guiding principle is taken, but also QCD-based and nonequilibrium-statistical arguments are considered. The rapid local equilibration of gluons and quarks in the initial stages of a relativistic heavy-ion collision is being modeled through exact analytical solutions of a nonlinear diffusion equation. On a similar time scale, stopping is accounted for in a QCD-inspired model, which is also incorporated in a nonequilibrium-statistical approach to compute the time evolution from the initial to the measured distribution functions of net protons. The production of charged mesons such as pions, kaons, and antiprotons is discussed in a phenomenological three-source relativistic diffusion model that emphasises the importance of the fragmentation distributions in addition to the usual fireball source. These are also shown to be essential in the phenomenon of limiting fragmentation, which had been confirmed experimentally at RHIC energies, and turns out to be consistent with the present LHC energies. The investigation of the dissociation of quarkonia in the QGP provides insights into the QGP properties, including an indirect measurement of its initial central temperature before the anisotropic expansion sets in. In asymmetric systems such as p-Pb at the current maximum LHC energy, the interplay of cold-matter and hot-medium effects has been studied, achieving a detailed understanding of the available data from the Large Hadron Collider. The more precise measurements from the forthcoming run 3 at the LHC are expected to provide deeper insights.

ACKNOWLEDGMENTS
I am grateful to members of UHD's Multiparticle Dynamics Group for discussions and collaborations, which resulted in common publications that are quoted in the references. Special thanks go to the Heidelberg students Viet Hung Dinh (now Orsay), Johannes Hölck, Benjamin Kellers, Niklas Rasch, Philipp Schulz, and Alessandro Simon.