Centrality dependence of limiting fragmentation

We investigate the centrality-dependent validity of the limiting-fragmentation hypothesis in relativistic heavy-ion collisions at energies reached at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC). A phenomenological analysis of Au-Au and Pb-Pb collisions within a three-source relativistic diffusion model (RDM) is performed at sqrt(s_NN) = 19.6, 62.4, 130, 200, 2760 and 5023 GeV using four centrality cuts at each energy. Linear and nonlinear expressions for the rapidity drift function are tested. Our results are compatible with the limiting-fragmentation conjecture for the investigated centralities in the full energy range. The number of particles in the fragmentation and fireball sources are found to depend on sqrt(s_NN) logarithmically and cubic-logarithmically, respectively.


Introduction
The occurence of limiting fragmentation (LF), or extended longitudinal scaling, had been predicted for hadron-hadron and electron-proton collisions by Benecke et al. [1]. It was first shown to be present in charged-hadron production at large pseudorapidities in the fragmentation region of pp data, in an energy range of √ s = 53 − 900 GeV [2]: The charged-particle pseudorapidity yield dN/dη does not depend on energy over a large range of pseudorapidities η = η − y beam , with the beam rapidity y beam . 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 characteristic feature of the particle production process, which turns out to be especially outstanding in relativistic heavy-ion collisions.
At the Relativistic Heavy Ion Collider RHIC in Brookhaven, limiting fragmentation was shown to occur in Au-Au collisions in the energy range √ s N N = 19.6 GeV to 200 GeV [3,4,5]. For a given centrality, the pseudorapidity distributions of produced charged particles were found to scale with energy according to the LF hypothesis in a given centrality class.
It is presently an open question whether limiting fragmentation will persist at the much higher incident energies of √ s N N = 2.76 and 5.02 TeV that are available at the CERN Large Hadron Collider (LHC) in Pb-Pb collisions since experimental results in the fragmentation region are not available due to the lack of a dedicated forward spectrometer. Nevertheless, it is most interesting to account for the collision dynamics more completely in this region. In Ref. [6], we had therefore studied central Pb-Pb collisions at LHC energies in a phenomenological model, with the result that LF can be expected to hold in central events.
We now extend this work to investigate the centrality dependence of limiting fragmentation in a three-source relativistic diffusion model (RDM, [7]). We consider equivalent centrality classes at four RHIC energies and two LHC energies in the range 0 − 30 %. At RHIC energies, a detailed comparison with PHOBOS data [8] is possible in all four centrality bins, whereas at LHC energies our analysis remains a model-dependent prediction. For all six energies and four centralities, we also deduce the number of produced charged hadrons in the respective two fragmentation sources and the fireball particle-production source, and determine their dependencies on √ s N N .
Our analysis complements microscopic approaches such as the multiphase transport model AMPT by Ko et al. [9] or HIJING [10,11] in order to assess whether centralitydependent LF is valid from RHIC to LHC energies. AMPT [9] had been tuned for the most central bin at RHIC and LHC energies [12]. In spite of disagreements with the LHC data in the midrapidity region, it had been concluded [13] 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 [14].
There exist also other phenomenological approaches such as the thermal model [15,16,17,18] and hydrodynamical models [19] which offer predictions regarding LF. The thermal model is appropriate to calculate particle production rates near midrapidity, but cannot be expected to predict distribution functions at forward rapidities that are needed to check the LF conjecture. It has nevertheless been employed in the forward region [20], concluding that LF should be violated at LHC energies. The relativistic diffusion model uses three sources for particle production: a midrapidity fireball source and two fragmentation sources [21,22,23,24,25,26]. The time evolution of the distribution functions is accounted for through solutions of a Fokker-Planck equation (FPE) for the rapidity variable which are subsequently transformed to pseudorapidity space through the appropriate Jacobian. The three sources can then be added to obtain the chargedhadron distribution that is used to check the LF hypothesis.
In Ref. [27] we have shown that a FPE in rapidity space can be derived from a nonequilibrium-statistical theory of non-Markovian processes in spacetime that are equivalent [28,29] to relativistic Markov processes in phase-space. The fluctuating background that is required for such a description to be valid is provided by the quarks and gluons in the fragmentation sources and in the fireball. A thermally equilibrated heat bath as in the theory of Brownian motion is not needed in the derivation. One obtains a FPE for time-dependent particle transport in rapidity space. Drift and diffusion terms are related through a fluctuationdissipation relation (FDR). With a constant diffusion coefficient and the FDR, the drift function in stopping can be calculated from the condition that the stationary solution of the FPE equals a distribution function derived in the color-glass condensate framework. For particle production, a similar path will be pursued.
In the present work about charged-hadron production and limiting fragmentation, we also use a FPE with constant diffusion coefficients, and either a linear dependence of the drift on rapidity y as in the original phenomenological RDM, or a sinh(y) dependence [30,6] that asymptotically leads to the Maxwell-Jüttner equilibrium distribution. In case of the nonlinear drift, the transport equation itself still remains linear, such that the time dependence of the three sources can be treated individually, and the corresponding results be added incoherently.
We briefly summarise the basic formulation in the next section. In sect. 3, we apply the model with sinh-drift term that requires a numerical solution to calculate chargedhadron pseudorapidity distributions in Au-Au and Pb-Pb at RHIC and LHC energies in four different centrality bins. In each case, the LF conjecture is tested. For central collisions, results are compared with analytical solutions obtained earlier for linear drift. In sect. 4 we determine the number of produced charged hadrons in the fragmentation sources and the fireball source at all six energies and four centralities, and discuss their energy dependence. The conclusions are drawn in sect. 5.

A phenomenological three-source model
The Lorentz-invariant cross section for produced particles in relativistic heavy-ion collisions is with the energy E = m T cosh(y), the transverse momentum p T = p 2 x + p 2 y , the transverse mass m T = m 2 + p 2 T , and the rapidity y = 0.5 ln[(E + p || )/(E − p || )].
In a three-source model for particle production, the rapidity distributions for all three sources k = 1, 2, 3 are obtained by integrating over the transverse mass m T where the normalisation constants c b k for the three sources depend on centrality, or impact parameter b. (Here and subsequently we omit the index b in all other variables such as N k ). The experimentally observable distribution dN/dy is evaluated in every centrality bin at the freezeout time, t = τ f , which corresponds to the interaction time τ int of Refs. [21,30]: 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 produced charged hadrons N k ch and adding them according to where 3 ≡ gg indicates that the fireball source is mostly arising from low-x gluon-gluon collisions.
Within the model with sinh-drift that we consider in this work, the superposition of particles from the three sources is still possible because the FPE is a linear partial differential equation. For symmetric systems, the problem is simplified by only considering the solution for the positive rapidity region and mirroring the result at y < 0. The parameters of the three-source model will be determined via χ 2 -minimisation with respect to the available data, and can be used in extrapolations and predictions.
In view of the high temperatures reached in relativistic collisions at RHIC and LHC energies, we rely on Boltzmann-Gibbs statistics and adopt the Maxwell-Jüttner distribution as the thermodynamic equilibrium distribution for In the relativistic diffusion model [21,22,25,27,30], the partial distribution functions R k (y, t) (k=1,2,gg) evolve in time towards this thermodynamic equilibrium distribution through solutions of the Fokker-Plannck equation with drift functions J k (y, t) and diffusion functions D k (y, t). The latter are taken to be constant coefficients D k in this work. If, in addition, the drift functions are assumed to 76 TeV Pb-Pb ALICE data [31] for charged-hadron production as in Refs. [7,6]. The solid curve is the overall distribution, the incoherent sum of the three sources. The distributions resulting from the fragments are symmetric (dotted). The fireball source is the essential contribution to the charged-hadron yield (dashed). (b) Numerical solution of the RDM with nonlinear drift. The distributions resulting from the fragments cover the full pseudorapidity range (dotted). The midrapidity source is wider and more pronounced compared to the linear-drift case (dashed).
be linearly dependent on the rapidity variable y, the FPE has the Ornstein-Uhlenbeck form [32] and can be solved analytically [21]. 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.
For constant diffusion and linear drift, the equilibrium limit of the FPE solution deviates, however, slightly from the Maxwell-Jüttner distribution, although the discrepancies are small and become visible only for sufficiently large times. In order to attain the correct stationary solution, the drift term must be modified according to [33,30] with a drift amplitude The drift force in the fragmentation sources k = 1, 2 grows with increasing distance in y-space from the beam rapidity, which enters through the initial conditions. The rapidity distribution at thermal equilibrium can then be derived [30] using Eqs. (2) and (4) as with C b being proportional to the overall number of produced charged hadrons N tot ch in the respective centrality bin. The actual distribution functions remain far from thermal equilibrium and the total particle number is evaluated based on the nonequilibrium solutions of the FPE, which are adjusted to the data in χ 2 -minimisations.
We determine the drift amplitudes A k in each centrality bin from the position of the fragmentation peaks. Diffusion coefficients can then be calculated as D k = A k T /m T from eq. (7). These refer, however, only to the diffusive processes. Since the fireball source and both fragmentation sources also expand collectively, the actual distribution functions are much broader [7] than what is obtained from eq. (7). We therefore use values for the diffusion coefficients (or the widths of the partial distributions) that are adapted to the data. The total particle number is then obtained in each centrality bin from the integral of the overall distribution function. We shall investigate four centrality bins at both, RHIC and LHC energies.
The RDM with linear drift has analytical solutions that can be used directly in χ 2 -minimisations with respect to the data, but numerical solutions of the FPE are required for the sinh-drift, as outlined in Ref. [30]. To arrive at a usable form for the computer, we transform the equation for R(y, t) into a dimensionless equation for f (y, τ ) through the introduction of a timescale t c that defines the dimensionless time variable τ = t/t c , such that To recover the drift and diffusion coefficients, one therefore has to specify a time scale. Since the drift determines the peak position, we choose the time-like variable τ such that the peak position of the experimental data is reproduced, leaving the diffusion strength γ as free parameter. Hence, for the three partial distributions in each centrality bin, there are three free parameters γ k , with the two values γ ( 1, 2) for the fragmentation sources being identical for symmetric systems such as Pb-Pb, but differing for asymmetric systems like p-Pb. The numerical solution is obtained using matlab's integration routine pdepe for solving parabolic-elliptic partial differential equations. We had shown in Ref. [30] that this method is very accurate when compared to results of finite-element methods such as DUNE [34] and FEniCS [35].

Charged-hadron production and limiting fragmentation
We insert physical values for the temperature T , the transverse mass m T , and the initial conditions f i (y, t = 0) in order to compare the model results to centrality-dependent data. Two distributions centered at the beam rapidities y beam = ± ln( √ s N N /m p ) 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 [30], we use gaussians with σ = 0.1. The same standard deviation is assumed 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, the critical value T = T cr = 160 MeV of the cross-over transition between hadronic matter and quark-gluon plasma is adopted. Experimental values are deduced for the transverse mass from measured transverse-momentum distributions.
We solve eq. (10) numerically for each centrality class and transform the results to rapidity distributions as discussed in Ref. [6] according to The constant C b is adjusted to the total number of produced charged hadrons in a given centrality bin b. At LHC energies, the fireball source yields the largest contribution to charged-hadron production. Particles that are produced from the fragmentation sources are not distinguishable from those originating from the fireball, but must be included in a phenomenological model. In particular, when regarding the limiting-fragmentation conjecture, the role of the fragmentation distributions is decisive since they determine the behavior of the distribution functions at large values of rapidity.
In case of unidentified charged particles, we first have to transform from rapidity-to pseudorapidity space in order to directly compare to data. The scattering angle θ determines the pseudorapidity variable η according to and we obtain the pseudorapidity distribution function dN dη from the rapidity distribution dN dy as dN dη = dy dη with the Jacobian for produced particles with mass m and 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. The full p T -distributions are, however, not available for all particle species that are included in the pseudorapidity measurements and hence, one has to make estimates. We had determined in Ref. [37] 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. Values at the other energies and centralities are found to vary between J 0 = 0.830 and 0.861. With eq. (14) for p T ≡ p eff T one obtains The mean mass m is calculated from the abundancies of pions, protons and kaons. Using J 0 , the Jacobian can be written independently from the values of m and p eff T as resulting in J (η) = cosh(η)[0.365 + cosh 2 (η)] −1/2 for central 2.76 TeV Pb-Pb collisions, and the same value at 5.02 TeV.
The effect of the Jacobian is most pronounced near midrapidity, where it generates the dip in the pseudorapidity distributions, as can be seen in fig. 1. Here, calculations in the RDM using both, linear drift [26] and sinh-drift [6] are compared with ALICE data for central Pb-Pb at 2.76 TeV [31]. The parameters and χ 2 -values for the sinh-drift are included in table 1, the linear-drift calculation is as in Ref. [26].
In this work, we are emphasizing the fragmentation region, where the Jacobian has almost no effect. We solve eq. (10) with sinh-drift using Dirichlet boundary conditions in four centrality classes at six incident energies, with parameters given in table 1. We perform χ 2 -minimizations with respect to the data in every centrality bin using Matlab. The resulting charged-hadron pseudorapidity distributions are shown in fig. 2 [31,36]. Due to the sinh-drift, the fragmentation distributions are less confined to the fragmentation region as in the linear-drift model, but extend over the whole pseudorapidity range. Hence, the Jacobian deforms also the fragmentation distributions in the midrapidity region, as shown already in fig. 1 for central collisions. Limiting fragmentation is clearly fulfilled within the RDM with sinh-drift in all four centrality classes when comparing with the Au-Au data at RHIC energies, and the results are consistent with LF for Pb-Pb data at LHC energies 2.76 and 5.02 TeV in all four centrality classes as well.

Charged-hadron content of the sources
As proposed in Ref. [7] for the RDM with linear drift, we now investigate the energy dependence of charged-hadron production in the three sources using the model with sinhdrift for each centrality class. Since the distributions resulting from the fragmentation sources have a different shape in the nonlinear-drift model, we can expect different results from the linear-drift case. The total number of produced charged hadrons follows a power law [7] The variable s is a dimensionless squared energy ratio s ≡ s N N /s 0 , with suitably chosen s 0 . The particle content produced by the fragmentation sources, N 1,2 depends logarithmically on s, The midrapidity source N gg , however, behaves differently [7]. Its width Γ is proportional to the beam rapidity y beam , Hence, the width scales logarithmically with s. The particles produced in the midrapidity source result mainly from low-x gluon-gluon interactions, with x the partonic longitudinal momentum fraction. The predicted cross section for such events is proportional to ln 2 s [38], satisfying the Froissart limit [39]. Since the cross section is directly proportional to the yield density, the midrapidity source distribution scales with ln s in width and with ln 2 s in yield density, such that the total functional dependence of the produced charged particles in the central source is expected to be [7,40] N gg ∝ ln 3 s.
We extract the values of N 1,2,gg from our RDM-analysis with nonlinear drift using suitable fit functions for the total number of produced charged particles, N tot = i N i , the number of particles produced in the fragmentation sources, N f = N 1 + N 2 , and the number of particles produced in the midrapidity source, N gg , as functions of the squared center of mass energy s N N . The fits are done in Matlab using the curve fit-function from the scipy.optimize package. The routine computes the best results via a χ 2minimization and also computes the covariance matrix, such that the square roots of its diagonal entries are the standard deviations. For the total number of charged particles a power law is expected. To be able to compare our findings directly with the results for linear drift in Ref. [7], we choose s 0 = 1 TeV 2 , such that s = s N N /1 TeV 2 and For the most central class (see table 2), we find a = (10.6 ± 0.2) × 10 3 , b = 0.229 ± 0.004 , consistent with the values a = 1.1 × 10 4 and b = 0.23 of Ref. [7] for linear drift. This is as expected, because the influence of the drift term on the total number of particles is negligibly small. The produced charged particles resulting from both fragmentation sources N f are with the results N 1,2 from section 3 for the individual sources. For the fragmentation sources we expect a logarithmic dependence Here we use a and s 0 as parameters. Since ln(s N N /s 0 ) = ln(s N N )−ln(s 0 ), s 0 determines the offset and should therefore be evaluated by the routine and not be fixed. This yields a = 313 ± 17, for the most central class. The values differ from the ones in the RDM with linear drift [7], because the latter tends to overestimate the effect of the fragmentation sources.
In the model with sinh-drift, the diffusion is so strong, that the fragmentation distributions spread over the whole pseudorapidity range. Hence, the form of the total distribution results in the case of a sinh-drift mainly from the central source -but the fragmentation sources are still relevant to decide whether LF is fulfilled. For the midrapidity source, the functional dependence on s is a cubic logarithm as discussed above, The results in the most central class are a = 9.1 ± 0.3, In Ref. [7] the RDM with linear drift resulted in a = 7.5, s 0 = 169 GeV 2 corresponding to a smaller yield in the central source. The difference is due to the overestimate of the fragmentation sources in the model with linear drift, which causes an underestimate of the midrapidity source. The results for the particle numbers in the centrality classes that we have investigated in the model with sinhdrift are shown in fig. 3. In central collisions, the results for the fragmentation sources using the model with linear drift as in Ref. [7] are also displayed, dot-dashed curve in (d). For the other centralities, analyses in the model with linear drift are not available. The midrapidity source becomes the main source for particle production at energies beyond the highest RHIC energy of 200 GeV, although at LHC energies the fragmentation sources still contribute substantially to charged-particle production.

Conclusion
We have investigated the centrality dependence of the limiting-fragmentation conjecture in charged-hadron pseudorapidity distributions in Au-Au and Pb-Pb collisions at RHIC and LHC energies. A three-source relativistic diffusion model with sinh-drift, which ensures the correct Maxwell-Jüttner equilibrium distribution, is the main basis of this study. It requires numerical solutions of the transport equation. For central collisions, a linear drift that allows for analytical solutions has also been tested.
The numbers of produced charged hadrons in the three sources as functions of incident energy have been deduced. The midrapidity source displays a cubic-logarithmic dependence on √ s N N and becomes the main source of particle production at LHC energies, but the fragmentation sources remain relevant and are essential to maintain limiting fragmentation. Our analysis shows that the RDM with three sources displays limiting-fragmentation scaling in agreement with the data in four centrality bins at RHIC energies. According to our results, limiting fragmentation is likely to be valid in the corresponding centrality bins at LHC energies, thus spanning a factor of almost 260 in collision energy.
This conclusion disagrees with expectations from simple parametrizations of the rapidity distributions, and also with predictions from the thermal model, which does not consider the fragmentation sources. The latter play an essential role in our nonequilibrium-statistical approach. It would be desirable that future upgrades of the LHC detectors will make it possible to test the limiting-fragmentation conjecture experimentally at LHC energies in different centrality bins.