Limiting fragmentation at LHC energies

We investigate the validity of the limiting-fragmentation hypothesis in relativistic heavy-ion collisions at energies reached at the Large Hadron Collider (LHC). A phenomenological analysis of central AuAu and PbPb collisions based on a three-source relativistic diffusion model (RDM) is used to extrapolate pseudorapidity distributions of produced charged hadrons from RHIC to LHC energies into the fragmentation region. Data in this region are not yet available at LHC energies, but our results are compatible with the limiting-fragmentation conjecture in the full energy range sqrt(s_NN) = 19.6 GeV to 5.02 TeV.


Introduction
The significance of the fragmentation region in relativistic heavy-ion collisions was realized when data on AuAu collisions in the energy range √ s NN = 19.6 GeV to 200 GeV became available at the Brookhaven Relativistic Heavy Ion Collider (RHIC) [1,2,3]. For a given centrality, the pseudorapidity distributions of produced charged particles were found to scale with energy according to the limiting fragmentation (LF), or extended longitudinal scaling, hypothesis: The charged-particle pseudorapidity distribution is energy independent over a large range of pseudorapidities η ′ = η − y beam , with the beam rapidity y beam . The existence of the phenomenon had been predicted for hadron-hadron and electron-proton collisions by Benecke et al. [4], and it was first shown to be present in pp data, in a range from 53 up to 900 GeV [5]. The fragmentation region grows in pseudorapidity with increasing collision energy and can cover more than half of the pseudorapidity range over which particle production occurs. The approach to a universal limiting curve is a remarkable feature of the particle production process, especially in relativistic heavy-ion collisions.
Presently it is not clear, however, whether limiting fragmentation will persist at the much higher incident energies that are available at the CERN Large Hadron Collider (LHC), namely, √ s NN = 2.76 and 5.02 TeV in PbPb collisions. Although detailed and precise ALICE data for charged-hadron production at various centralities are available in the midrapidity region [6,7] for both incident energies, experimental results in the fragmentation region are not available due to the lack of a dedicated forward spectrometer. This region is, however, most interesting if one wants to account for the collision dynamics more completely. In this work, we investigate to what extent limiting fragmentation can be expected to occur in heavy-ion collisions at LHC energies. Given the lack of LHC data in the fragmentation regions, one has to rely on either microscopic approaches such as the multiphase transport model AMPT by Ko et al. [8] or HIJING [9,10] in order to assess whether LF is valid at LHC energies, or on phenomenological models. Among these are the thermal model [11,12,13,14], hydrodynamical approaches [15], or the relativistic diffusion model (RDM) with three sources for particle production: a midrapidity source and two fragmentation sources [16,17,18,19,20,21]. Here, the time evolution of the distribution functions is accounted for through solutions of a Fokker-Planck equation (FPE) for the rapidity variable which is subsequently transformed to pseudorapidity space through the appropriate Jacobian.
Regarding microscopic approaches, the AMPT code [8] had been tuned for the most central bin at RHIC and LHC energies in Ref. [22]. There are some disagreements with the LHC data in the midrapidity region, but AMPT is in accord with longitudinal scaling at LHC energies. In Ref. [23] it had already been concluded that AMPT and other microscopic codes reproduce LF at RHIC energies. The same conclusion had been drawn from calculations in the color-glass-condensate framework [24].
The ALICE collaboration has argued in Ref. [25], in accordance with our results in Ref. [26], that their 2.76 TeV PbPb data are in agreement with the validity of extended longitudinal scaling -within the uncertainties which arise mainly from the extrapolation of the charged-particle pseudorapidity density from the measured region to the rapidity region of the projectile where no data are available. In their analysis, the extrapolation into the forward η-region was done using the difference of two gaussian functions as detailed in Ref. [25]. Both Gaussians are centered at midrapidity, and the second (subtracted) Gaussian simulates the central dip that is mostly due to the jacobian transformation from y-to η-space.
This procedure leads ALICE to conclude that the 2.76 TeV PbPb data are consistent with the LF hypothesis. Related, but different, extrapolation schemes give similar conclusions. As an example, we have fitted the ALICE midrapidity data with a sum of two Gaussians that are peaked at the experimental maxima, used the proper jacobian transformation from yto η-space, and determined the corresponding parameters for PbPb collisions in a χ 2 -minimization. Again, the resulting pseudorapidity distribution functions fulfill the LF hypothesis at both LHC energies, 2.76 and 5.02 TeV. However, such an extrapolation procedure is a rather arbitrary scheme without any physical basis. In contrast, the thermal model [11,12,13,14] has a macroscopic physical basis, which is appropriate to predict particle production rates at midrapidity -but it is questionable whether it is suited to predict distribution functions, in particular at forward rapidity. Still, it has been used in Ref. [27] in the forward region, with the conclusion that limiting fragmentation should be violated at LHC energies.
In this work, we investigate whether limiting fragmentation in heavyion collisions at LHC energies can be expected to be fulfilled in yet another phenomenological model, the three-source relativistic diffusion model (RDM) [16,20,21]. We briefly summarize the basic formulation in the next section. In section 3, we apply the model both in its analytically solvable version with linear drift, and with a sinh-drift term that requires a numerical solution, to calculate net-proton rapidity distributions at RHIC energies. Here, only the fragmentation sources contribute, and thus can be directly compared to data. In section 4 we apply the model including a central source to charged-hadron production in central AuAu and PbPb collisions at RHIC and LHC energies, in order to test whether limiting fragmentation is fulfilled at LHC energies. The conclusions are drawn in section 5.

A phenomenological three-source model
In relativistic heavy-ion collisions, the relevant observable in stopping and particle production is the Lorentz-invariant cross section with the energy E = m ⊥ cosh(y), the transverse momentum p ⊥ = p 2 x + p 2 y , the transverse mass m ⊥ = m 2 + p 2 ⊥ , and the rapidity y. We shall first investigate in the next section rapidity distributions of protons minus produced antiprotons, which are indicative of the stopping process as described phenomenologically in a relativistic two-source diffusion model (RDM) [16,28] or in a QCD-based approach [29]. This motivates the relevance of the fragmentation sources not only in stopping, but also in particle production at relativistic energies. In section 4, we then switch to a threesource model for particle production, with the fireball source contributing most of the produced charged hadrons at RHIC and LHC energies when compared to each of the fragmentation sources. This central source does not contribute to stopping because particles and antiparticles are produced in equal amounts. The three-source model is visualized schematically in Fig. 1 for a symmetric system such as AuAu or PbPb.
The rapidity distributions for all three sources k = 1, 2, 3 are obtained by integrating over the transverse mass with normalization constants c k that depend on the number of participants at a given centrality for k = 1, 2. The experimentally observable distribution dN/dy is evaluated in the time-dependent model at the freeze-out Figure 1: (color online) Schematic representation of the three-source model for particle production in relativistic heavy-ion collisions at RHIC and LHC energies in the center-ofmass system: Following the collision 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.
time, t = τ f . The latter can be identified with the interaction time τ int of Refs. [16,28]: the time during which the system interacts strongly. The full rapidity distribution function for produced charged hadrons is obtained by weighting the three partial distribution functions with the respective numbers of particles, and adding them incoherently: where the index 3 ≡ gg is meant to emphasize that the fireball source is mostly arising microscopically from low-x gluon-gluon collisions. The incoherent addition of the three sources applies also to the model with sinh-drift that we consider in this work, because the FPE is a linear partial differential equation, allowing for linear superposition of independent solutions. For a symmetric system, one can further simplify the problem by only considering the solution for the positive rapidity region and mirroring the result at y < 0.
The parameters of the three-source model -which will be detailed in the following -are then determined via χ 2 -minimization with respect to the available data, and can be used in extrapolations and predictions [21]. In stopping, the relevant distribution function is given by the incoherent sum of the fragmentation sources only, We rely on Boltzmann-Gibbs statistics and hence adopt the Maxwell-Jüttner distribution as the thermodynamic equilibrium distribution for t → ∞ The nonequilibrium evolution of all three partial distribution functions R k (y, t) (k = 1, 2, gg) towards this thermodynamic equilibrium distribution is accounted for in the relativistic diffusion model [16,17,20,28] through solutions of the Fokker-Planck equation with suitably chosen drift functions J k (y, t) and diffusion functions D k (y, t).
If the latter is taken as a constant diffusion coefficient D k , and the drift function assumed to be linearly dependent on the rapidity variable y, the FPE has the Ornstein-Uhlenbeck form [30] and can be solved analytically [16]. For t → ∞ all three subdistributions approach a single Gaussian in rapidity space which is centered at midrapidity y = 0 for symmetric systems, or at the appropriate equilibrium value y = y eq for asymmetric systems.
In case of stopping, only the two fragmentation distributions contribute, approaching the thermal equilibrium distribution for t → ∞, as will be shown in the next section. The equilibrium limit of the FPE solution for constant diffusion and linear drift is, however, found to deviate slightly from the Maxwell-Jüttner distribution. Although the discrepancies are small and become visible only for sufficiently large times, we use the RDM with the sinh-drift which ensures that the solution for t → ∞ yields the Maxwell-Jüttner distribution Eq. (5), as was discussed in Refs. [31,32]. This induces a special form of the fluctuation-dissipation theorem (FDT) that connects diffusive and dissipative phenomena in the collision, namely With Eqs. (2) and (5), the rapidity distribution at thermal equilibrium can then be derived [28] as 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 χ 2minimizations. 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 ⊥ . These refer, however, only to the diffusive processes, and since the fireball source and both fragmentation sources also expand collectively, the actual distribution functions will be much broader than what is obtained from Eq. (8). Hence, we shall use values for the diffusion coefficients (or the widths of the partial distributions) that are adapted to the data in both stopping and particle production. The total particle number is then obtained from the integral of the overall distribution function.
Whereas the RDM with linear drift has analytical solutions that can be used directly in χ 2 -minimizations of the data, numerical solutions of the FPE are required for the sinh-drift, as described in Refs. [28,32]. To arrive at a usable form for the computer, we transform the equation for R(y, t) into its dimensionless version for f (y, τ ) by introducing a timescale t c , defining the dimensionless time variable τ = t/t c . It follows that ∂ ∂t = ∂ ∂τ t −1 c and hence ∂f ∂τ Since The result is the dimensionless Eq. (11) depending only on the ratio γ = T /m ⊥ of temperature T and transverse mass m ⊥ which is a measure of the strength of the diffusion, To recover the drift and diffusion coefficients, one has to specify a time scale (or the other way round). Considering that it is only the drift term that is responsible for determining the peak position, we choose the time-like variable τ such that the peak position of the experimental data is reproduced. This leaves the diffusion strength γ as free parameter. In case of three partial distributions, there are three free parameters γ k . Here, the two values for the fragmentation sources are identical for symmetric systems such as PbPb, but differ for asymmetric systems like pPb.
We calculate the numerical solution using matlab's integration routine pdepe for solving parabolic-elliptic partial differential equations. It was shown in Ref. [28] that this method is very accurate when compared to results of finite-element methods such as DUNE [33] and FEniCS [34].
To compare the simulation to experimental data, we have to insert relevant values for T , m ⊥ , and the initial conditions. The beam rapidity y beam is determined by the center-of-mass energy per nucleon pair as y beam = ± ln( √ s N N /m p ). Two gaussian distributions centered at the beam rapidities with a small width that corresponds to the Fermi motion represent the incoming ions before the collision. The exact width of the initial distribution does not have a large effect on the time evolution [28]; here we use σ = 0.1. The same standard deviation is taken for the initial condition of the midrapidity source, which is centered at y = 0 for a symmetric system, and at y = y eq for asymmetric systems.
For the temperature, we take the critical value T = T cr = 160 MeV for the cross-over transition between hadronic matter and quark-gluon plasma. Regarding the transverse mass, experimental values are deduced from measured transverse-momentum distributions.
The constant C is chosen in case of stopping such that the total number of particles corresponds to the number of participant protons in the respective centrality bin. In particle production, C is adjusted to the total number of produced particles for a given centrality.

Fragmentation sources in stopping
To emphasize the relevance of the fragmentation sources, we first investigate stopping and calculate net-proton rapidity distributions in central AuAu collisions at RHIC energies of 200 GeV, where data are available from Ref. [35]. Theoretical calculations are usually performed for net-baryon distributions [29] because the total baryon number is a conserved quantity, but since experimentally only net-proton distributions are available, we convert to net-protons via Z/A = 79/197 = 0.40. As discussed, the central source cancels out in stopping, because particles and antiparticles are produced in equal amounts. Rather precise net-proton results had been obtained at SPS energies for central PbPb at √ s NN = 17.3 GeV, where the NA49 collaboration succeeded to measure across the fragmentation peak in a fixed-target experiment [36]. These results can be well reproduced in the linear RDM [37], and also in a QCD-based approach [29], but here we are interested in higher energies, namely, the RHIC and LHC region. For AuAu at a RHIC energy of 200 GeV, it was not possible to measure the fragmentation peaks, because the forward spectrometer hit the beam pipe at large rapidities. At LHC energies, a forward spectrometer with particle identification in the region of the expected fragmentation peaks [29] is not available. Still, the BRAHMS stopping data of Ref. [35] shown in Fig. 2 indicate the rise towards the fragmentation peaks, which was later corroborated by more recentalbeit preliminary -data near the peak region [39]. We now compare these AuAu RHIC data with the fragmentation distributions that arise from the three-source model with both linear and sinhdrift. For AuAu the transverse mass m ⊥ is obtained from the p ⊥ -spectra using m ⊥ = m 2 p + p 2 ⊥ and Eq. (1) for the yields. As proposed in Ref. [35], Gaussians are fitted to the invariant yields [28]. The results for protons and antiprotons are averaged to obtain The results of the RDM-calculation with the sinh-drift are shown as dashed curves in Fig. 2. The dimensionless time parameter has been adjusted as τ = 0.08 with the above value of the diffusion strenght γ to yield a fragmentation-peak position of y = ±3.1, in accordance with the data of Refs. [35,39] (only the final data of Ref. [35] are shown here). The calculated distribution function is, however, by far too narrow, because the theoretical expression from Eq. (8) does not account for collective expansion. In Ref. [40], the longitudinal expansion velocity v || had actually been calculated from the difference between the theoretical distribution function, and the data. The solid curve represents a numerical solution of the FPE with adapted diffusion strength γ = 33, it clearly shows the fragmentation peaks. 6 GeV (lower data points, [38]) are consistent with the limiting fragmentation hypothesis (LF). The total density distribution can be separated into three parts, one resulting from the midrapidity source (dashed curve) and two from the fragmentation sources (dot-dashed curves).
In Fig. 2 we also display the corresponding equilibrium solutions, which are centered at midrapidity for symmetric systems. We use the theoretical FDT-value γ = 0.139, and display the numerical solution of the FPE with sinh-drift for τ → ∞ as dot-dashed curve. It agrees with the Maxwell-Jüttner distribution Eq. (5) for m ⊥ = 1.15 GeV and T = 160 MeV (crosses). For comparison, the equilibrium result of the RDM with linear drift is also shown (triangles). It is a Gaussian with a variance σ 2 eq = γ = t c D = A −1 D = T /m ⊥ = 0.139 corresponding to a width Γ FWHM = √ 8γ ln 2 = 0.88. The normalization C is such that the integral of the total distribution yields the number of participant protons in 0−5% central AuAu collisions, N p = N B Z/A = 357×79/197 ≃ 143, with the number of participant baryons N B = 357 ± 8 from a Glauber calculation [35]. The resulting distribution function deviates only slightly from the Maxwell-Jüttner equilibrium distribution.
Equilibrium distributions with physical values for the diffusion strengths that include collective expansion would be much broader, but they do not exhibit fragmentation peaks with a midrapidity valley, and hence, equilibrium models are not suited to describe stopping distributions.

Particle production and limiting fragmentation
In charged-hadron production, we consider the sum of produced charged particles and antiparticles. Hence, the fireball source has to be added, and yields the essential contribution to charged-hadron production in heavy-ion collisions at RHIC and LHC energies. Particles that are produced from the fragmentation sources are not directly distinguishable from those originating from the fireball, but still the fragmentation sources are relevant and must be included in a phenomenological model. In particular, when regarding the limiting-fragmentation conjecture, the role of the fragmentation distributions will turn out to be decisive since they determine the behavior of the rapidity distribution functions at large values of rapidity.
For unidentified charged particles, we first have to transform from rapidityto pseudorapidity space in order to directly compare to data. The pseudorapidity variable η is uniquely determined by the scattering angle θ and the pseudorapidity distribution function dN dη is obtained from the rapidity distribution dN dy through the transformation with the Jacobian for produced particles with mass m and transverse momentum p ⊥ . The transformation depends on the squared ratio (m/p ⊥ ) 2 of mass and transverse momentum of the produced particles. Hence, its effect increases with the mass of the particles and is most pronounced at small transverse momenta. In principle, one has to consider the full p ⊥ -distributions, which are, however, not available for all particle species that are included in the pseudorapidity measurements. In Ref. [26], we have determined the Jacobian J 0 at η = y = 0 in central 2.76 TeV PbPb collisions for identified π − , K − , and antiprotons from the experimental values dN dη | exp and dN dy | exp as J 0 = 0.856. Solving Eq. (17) for p ⊥ ≡ p eff ⊥ yields with a mean mass m which may be calculated from the abundancies of pions, protons and kaons. With the introduction of J 0 , the Jacobian can then be written independently from the values of m and p eff ⊥ as which results in J(η) = cosh(η)[1.365 + sinh 2 (η)] −1/2 for central 2.76 TeV PbPb collisions. The effect of the Jacobian is most pronounced near midrapidity, where it is essential to generate the dip in the pseudorapidity distributions, as is obvious from Fig. 3: A calculation in the RDM with linear drift [21] is compared with ALICE data for central PbPb at 2.76 TeV [6], with five parameters and χ 2 -values from Tab. 1. The optimization is done using Python. The Jacobian has almost no effect in the fragmentation region, which we are emphasizing in this work. In Fig. 3, we also compare the RDM-solution for 2.76 TeV PbPb with central AuAu data [38] at √ s NN = 19.6 GeV from the PHOBOS collaboration at RHIC. When plotted as function of η − y beam , limiting fragmentation is obviously fulfilled in the relativistic diffusion model with linear drift. This would not be the case if only the central fireball source was considered, as limiting fragmentation is a consequence of the appearance of the fragmentation sources. This result indicates that limiting fragmentation can be fulfilled from RHIC to low LHC energies in the relativistic diffusion model.  We now proceed to investigate the consequences of the model with sinhdrift, with emphasis on the fragmentation region. We solve Eq. (11) for central PbPb collisions at √ s NN = 2.76 and 5.02 TeV using Dirichlet boundary conditions, with parameters given in Tab. 2, and the same Jacobian for both energies. The resulting charged-hadron pseudorapidity distributions are shown in Fig. 4, again plotted as functions of η − y beam together with central AuAu data [38] at √ s NN = 19.6 GeV from the PHOBOS collaboration at RHIC. As a consequence of the sinh-drift, the fragmentation distributions are now much less confined to the fragmentation region, but extend into the whole pseudorapidity range that is accessible for produced charged hadrons. Hence, the Jacobian deforms also the fragmentation distributions in the midrapidity region. In the fragmentation region, limiting fragmentation is very well fulfilled when comparing the results at LHC energies to the 19.6 GeV AuAu data. Our overall results from the relativistic diffusion model with linear and sinh-drift are summarized in   [6,7]. The midrapidity sources (dashed curves) remain symmetric, but the fragmentation sources (dotted and dot-dashed curves) are asymmetric due to the sinh-drift. Comparison with the 19.6 GeV AuAu RHIC-data [38] confirms the consistency with the limiting fragmentation hypothesis in both cases.  [6,7], and PHOBOS data [38]. The zoom into this region shows that the RDM with sinh-drift is consistent with limiting fragmentation at RHIC and LHC energies.
200 GeV [38] with ALICE data [6,7] and our results for central PbPb at LHC energies of 2.76 and 5.02 TeV from the relativistic diffusion model with both linear and sinh-drift. The parameters are summarized in Tab. 1 and Tab. 2. As expected, the RDM with sinh-drift (solid curves) gives a better representation of the data in the fragmentation region as compared to the analytical linear model (dashed and dot-dashed curves). The inset shows that the RDM with sinh-drift is in agreement with the limiting-fragmentation conjecture at the available RHIC and LHC energies.

Conclusion
Our analysis of charged-hadron pseudorapidity distributions in central PbPb collisions within a three-source relativistic diffusion model indicates that the phenomenon of limiting fragmentation scaling can be expected to hold at RHIC and LHC energies, spanning a factor of almost 260 in collision energy. This conclusion is in line with results from microscopic numerical models such as AMPT, but it disagrees with expectations from simple parametrizations of the rapidity distributions such as the difference of two Gaussians, and also with predictions from the thermal model. The latter does not explicitly treat the fragmentation sources, the model refers only to particles produced from the hot fireball. In our phenomenological approach that confirms the validity of the LF hypothesis at LHC energies, however, the fragmentation sources play an essential role. It remains to be seen whether future upgrades of the detectors will make it possible to actually test the limiting-fragmentation conjecture experimentally at LHC energies.