Testing nonextensive statistics in relativistic heavy-ion collisions

Numerical solutions of the nonlinear Fokker-Planck equation (FPE) which has been associated with nonextensive q-statistics show that the available data on rapidity distributions for stopping in relativistic heavy-ion collisions cannot be reproduced with any permitted value of the nonextensivity parameter (1 < q < 1.5). This casts doubt on the nonextensivity concept that is widely used in relativistic heavy-ion physics.


Introduction
Nonextensive statistics proposes an extension of Boltzmann statistics through the concept of a non-additive q-entropy. It has been used in a nonlinear Fokker-Planck equation (FPE) for rapidity distributions, and applied to calculate rapidity and transverse momentum distributions for produced and stopped charged particles in relativistic heavy-ion collisions. In the corresponding experiments, the measured charged-hadron rapidity distributions are found to be very broad compared to thermal model predictions [1], and the discrepancy increases strongly with energy. This finding, as well as correspondingly broad net-proton (proton minus antiproton, or stopping) distributions [2,3], indicates thermal diffusion plus collective expansion.
Both effects may be accounted for phenomenologically in a linear diffusion model [4] with expansion, or else in abundant hydrodynamic approaches (see e.g. [5] for a review). It has, however, been stipulated that the so-called nonextensive q-statistics as proposed by Tsallis et al. [6] can simultaneously account for thermal and collective effects merely through a suitable choice of q [7][8][9].
With values of 1 < q < 1.5, the Fokker-Planck equation that has been used to model rapidity distributions [4] becomes nonlinear, it has an exponent (2 − q) in the diffusion term [6][7][8][9]. This is supposed to account for long-range forces that cause collective expansion, and is considered to be a fundamental property of the system like the temperature T . It goes along with a modified definition of the system's entropy [6] which is, however, controversial on fundamental grounds [10,11]. This approach is quoted in [7] as having the additional reward that the Einstein relation between drift and diffusion coefficient -that is valid in the theory of Brownian motion -could be maintained.
Indeed it has been claimed in [7,8] that such a procedure provides fits to stopping data in PbPb and AuAu collisions at energies reached at the Super Proton Synchrotron (SPS) and the Relativistic Heavy Ion Collider (RHIC) with values 1 < q < 1.5. (No stopping distributions will be available in the foreseeable future from the Large Hadron Collider (LHC) due to the lack of a suitable forward spectrometer). However, the solutions of the nonlinear FPE in these calculations are obtained by starting from the stationary solution as proposed in [6], and then solving the problem for time-dependent temperatures T (t) and mean values of the rapidity y m (t).
It is the aim of this work to solve the nonlinear FPE directly with realistic physical initial conditions. We shall use several independent numerical schemes, but without a predetermined form of the solutions -such as taking the form of the stationary solution as a basis for the time-dependent case as had been done in [7] -, and then try to fit the measured stopping distributions at SPS and RHIC energies with a value of q > 1.
Some relevant model ingredients for the linear and nonlinear cases are summarized in the next section. The numerical calculations are prepared and tested in the subsequent section, and compared with stopping data for PbPb and AuAu collisions at SPS and RHIC energies. We show that it is not possible to fit the data using solutions of the nonlinear FPE with values of the nonextensivity coefficient 1 < q < 1.5. Instead we fit the data in the linear model [12] with an adjusted diffusion coefficient to account for both, nonequilibrium thermal broadening, and collective expansion. The results are then briefly summarized. Preceding this proceedings article, they have been published in Ref. [13].

Basic considerations
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. In this work, we concentrate on rapidity distributions of protons minus produced antiprotons, which are indicative of the stopping process as described phenomenologically in a relativistic diffusion model (RDM) [4,12], or in a QCD-based approach [14]. The rapidity distribution is then obtained by integrating over the transverse mass dN dy with a normalization constant C that depends on the number of participants at a given centrality. The experimentally observable distribution dN/dy is calculated for the freeze-out time, t = τ f . The latter can be identified with the interaction time t = τ int of Refs. [4,12]: the time during which the system interacts strongly. We rely on Boltzmann-Gibbs statistics and hence, adopt the Maxwell-Jüttner distribution as the thermodynamic equilibrium distribution for t → ∞ The Boltzmann-Gibbs definition of the entropy is S = −k B Ω i=1 p i ln(p i ), where p i equals the probability of the system to be in the microstate i. In the case of equal probabilities and a total number of states Ω it follows that p i = p = 1 Ω , and (with k B ≡ 1) which is the well-known expression for the entropy. It is an extensive quantity: It becomes larger as the size of the system increases. To show this, take takes two systems A and B which do not interact. The number of available microstates in the combined system is equal to the product of the ones in the individual systems as they do not interact, Ω(A + B) = Ω(A) Ω(B).
Inserting this into the definition of entropy, one gets Although classical thermodynamics is a very successful theory, discrepancies with respect to data can arise. This is particularly relevant in the case of nonequilibrium systems, such as relativistic heavy-ion collisions. For such systems, various different concepts of entropy beyond Boltzmann-Gibbs have also been developed. In particular, Tsallis has proposed to resort to nonextensive statistics [6,15] where the entropy does not fulfill Eq. (4) but is instead given by with the entropic index q ∈ R. Here, the logarithm which causes the additivity of the entropy has been replaced by the non-additive q-logarithm , and q measures the degree of nonextensivity. The inverse of the q-logarithm is the q-exponential e x q that solves the differential equation dy/dx = y q through In the limit q → 1, S q is equal to S because provided the last term in Eq. (7) is neglected, There is, however, no clearly defined physical process that would warrant a generalization from S to S q , and no theory available to calculate the nonextensivity exponent q from first principles. It can still successfully be used as an additional fit parameter, in particular for p ⊥ -distributions in pp and AA collisions at relativistic energies which show a transition from exponential to power-law behaviour that the e x q -function properly describes with q ∈ (1, 1.5). From a more fundamental point of view, the approach is controversial [10,11]. In this work, we test its applicability to rapidity distributions in relativistic heavy-ion collisions.

Fokker-Planck equation
The general form of the linear Fokker-Planck equation (FPE) is [16] where J is called the drift coefficient and D is the diffusion coefficient. Here we denote the independent variable as y because it will later considered to be the rapidity. Even for coefficients J and D that are not time dependent it is generally difficult -if not impossible -to find analytical solutions. Two important analytically solvable examples are J(y, t) = 0, D(y, t) = D (Wiener process) and J(y, t) = −αy, D(y, t) = D (Uhlenbeck-Ornstein process [17]). For more complicated problems, numerical methods are employed.
In the relativistic diffusion model, the time evolution of the rapidity spectra has been modeled by a FPE. At first an Uhlenbeck-Ornstein (UO) ansatz has been tested in Ref. [4]. The stationary solution in such a case is determined as C has to be equal to zero because otherwise W < 0, This does not correspond to the equilibrium distribution from Eq. (3) and therefore another drift term is needed. We see from the above calculation that a stationary solution W ∝ e −V(y) results from a drift term V (y). With the drift J(y) = −A sinh(y) one gets the desired stationary solution [7,12] with A = m ⊥ D T , which can be interpreted as a fluctuation-dissipation relation similar to one known from Brownian motion, D = bT with the mobility b. Hence, the dissipation as described by the amplitude of the drift term can be related to the diffusion coefficient that is responsible for the fluctuations.
This particular sinh-drift term has also been investigated in Ref. [12] and the result was -as in the simple UO model [18] -that the fluctuation-dissipation relation is violated: The diffusion is too small to account for the experimental data.
The canonical interpretation of this result is that collective expansion occurs in the quarkgluon-plasma phase and enhances the width. One way to match the observation is to increase the diffusion coefficient, attributing the effect to collective expansion [12,18]. Indeed a general form of the fluctuation-dissipation theorem has been used in relativistic hydrodynamic calculations that describe systems exhibiting longitudinal collective expansion [19].
Within the Fokker-Planck framework, another possibility is to change the underlying equation in order to account for the 'anomalous' diffusion [7][8][9]. In the latter approach which we want to test here, one extends the model to a nonlinear FPE Analytical solution strategies for this equation in case of ν 1, µ 1 are not readily available. However, one can connect Eq. (12) with the nonextensive entropy Eq. (5). Indeed, Tsallis and Bukman have shown in [6] that the result of maximizing the entropic form leads to the function When assuming a drift term J(y, t) = −αy and a constant diffusion coefficient D(y, t) = D, the function p q (y, t) ≡ W(y, t) solves the partial differential equation Eq. (12) with additional conditions on β(t), y m (t) and Z q (t). One can identify q from the entropic form with the exponents µ and ν of Eq. (12) as q = 1 + µ − ν [6].
This identification is actually only justified in the case of the above linear drift, which is not the one we will use because the Boltzmann equilibrium form requires a sinh-drift. It was also shown in Ref. [6] that in order to conserve the norm, µ = 1 is required, and since we model a probability distribution we set µ to one, such that the exponent of the diffusion term becomes ν = 2 − q.

Numerical results and comparison with experimental data
To arrive at a usable form for the computer, we transform the equation for W(y, t) into its dimensionless version for f (y, t) by introducing a new timescale t c , resulting in the dimensionless time variable τ = t/t c . It follows that ∂ ∂t = ∂ ∂τ t −1 c and further Since The result is the dimensionless Eq. (15) 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 get physical values for 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 are free to chose the time τ such that the peak position of the experimental data is reproduced. This leaves as free parameters the diffusion strength γ and the nonextensivity parameter q. We calculate the solution using two different methods in order to gain insight about the accuracy. The more straightforward one was using matlab's integration routines for solving parabolic-elliptic PDEs. The second, more elaborate method, was implementing it in a finiteelement-method framework (FEM) (DUNE [20] and FEniCS [21]), see Ref. [13] for further details.
To compare the simulation to experimental data, we have to insert relevant values for T , m ⊥ and the initial conditions, most importantly y 0 . The value of the beam rapidity y 0 is determined by the center-of-mass energy per nucleon pair as y 0 = ln( √ s NN /m p ). Two Gaussian distributions centered at ±y 0 with a small width σ that corresponds to the Fermi motion represent the incoming ions before the collision. The exact value of σ does not have a large effect on the time evolution [12]; here we use a value of 0.1. For the temperature, we take the critical value 160 MeV for the transition between hadronic matter and quark-gluon plasma. The actual freeze-out temperature is smaller (T = 118 ± 5 MeV for PbPb at SPS energies [2]); overestimating the temperature will increase the diffusion. For 17.2 GeV PbPb, the transverse mass is taken to be m ⊥ = 1.17 GeV as the average transverse momentum p ⊥ is around 0.7 GeV [2]. The dimensionless diffusion strength γ is thus 0.137. Corresponding values for 200 GeV AuAu are given in Ref. [13].
The results are then transformed to a rapidity distribution [12]. Rewriting Eq. (2) and replacing d 3 N/dp 3 with the computed distribution f (y, t), we obtain The constant C is chosen such that the total number of particles for 0 − 5% corresponds to the number of participant protons in this centrality bin.
In order to check the numerical implementation, we compare it to analytically solvable problems. At first, we consider the UO model and compare the numerical solution of Eq. (15) for different values of q, first for q = 1 (where both should be the same) and then for other values of q. This gives a first idea about the impact of the non-linearity parameter on the evolution.
The numerical result for q = 1 is identical with the analytical solution, which validates the numerical method. By increasing q, the peaks are slightly smeared out, giving an overall flatter shape than before. This is expected since a larger diffusion coefficient will spread out the profile faster. As the next step, we consider the problem solved analytically by Borland et al. in Ref. [22] ∂ f ∂t The solution assumes that the initial condition is functionally equal to the stationary solution, except for time-dependent coefficients. In particular, both the stationary solution, and the initial conditions are centered at y = 0, which is essential to obtain the analytical solution of the time-dependent problem. The agreement between analytical and numerical solution (see Fig. 1) in the case of initial conditions that are centered at y = 0 further supports the correctness of the implementation.
In the case of a heavy-ion collision, however, the initial distributions are both off-center at the values of the beam rapidities, whereas the stationary solution that is obtained for t → ∞ is centered at rapidity y = 0 in symmetric systems. Hence, the Borland et al. analytical solution cannot be used: The time-dependent equation must be solved with the beam rapidities y beam = ±y 0 defining the initial conditions (δ-functions, or Gaussians with a width that is determined by the Fermi motion), and the solution in the heavy-ion case drifts with increasing time towards midrapidity. Although the Borland solution cannot describe our physical situation, it offers the possibility to compare the anomalous diffusion to an analytic solution. GeV: While a larger q does broaden the distribution, the effect is by far too small to come close to the experimental results. In order to reproduce the measured data for PbPb, we have to adopt a diffusion strength of around 1.5 while the one predicted by the fluctuationdissipation relation is around 0.137, the difference being a factor 11, see the upper frame of Fig. 3, for q = 1. As already mentioned, such a large enhancement in the required broadening cannot be compensated by the proposed nonlinearity due to q-statistics.
The comparison with AuAu stopping data at the maximum energy of 200 GeV reached at RHIC shows that here the discrepancy between the diffusion strength from the fluctuationdissipation relation (γ = 0.12) and the one required to fit the data (γ = 8) with an adjusted value of time is even larger, see lower frame of Fig. 3. This means that introducing a nonlinearity into the diffusion term cannot account for the observed rapidity spectra at SPS and RHIC energies. Since the widths are too narrow, there has to be an additional expansion process that takes place during the reaction that cannot be accounted for by q-statistics. This result is in obvious contrast to the findings of Refs. [7][8][9], where an approximate solution of Eq. (15) had been used. Rather than adopting nonextensive statistics, we therefore employ the model with linear diffusion q = 1, and the drift term imposed by the stationary solution [12]. We find physical values for the drift and diffusion coefficients in stopping using the two data sets from NA49 [2] at √ s NN = 17.2 GeV, and from BRAHMS [3] at √ s NN = 200 GeV.
The failure to interpret the broad rapidity distributions observed in the stopping process of relativistic heavy-ion collisions within q-statistics refers specifically to the solution of the nonlinear Fokker-Planck equation Eq. (12) which arises within nonextensive statistics [6,15]. Among the abundant publications that compare q-statistics with data in particular for p ⊥distributions (e.g. [23][24][25]), but also for rapidity distributions ( [24]), only few such as [7,8] refer to an explicit -but approximate -solution of the basic nonlinear FPE.

Conclusion
We have tested the nonextensive paradigm in a well-defined application to rapidity distributions for stopping in relativistic heavy-ion collisions. When using first Boltzmann-Gibbs statistics with a linear drift, or with the proper sinh-drift that secures the correct Maxwell-Jüttner equilibrium limit, the calculated widths of the rapidity distributions are by far too small compared to the data. The usual explanation for the discrepancy is that collective expansion must be taken into account that broadens the distributions far beyond the statistical widths.
Nonextensive statistics claims that it implicitly accounts for long-range interactions that cause expansion and hence, it should be feasible to reproduce the experimental widths through an adjustment of q > 1. In this context a nonlinear Fokker-Planck equation is available which had been solved previously by Lavagno et al. [7][8][9] using an approximate solution scheme. However, when actually solving this equation numerically, this turns out to be impossible, as shown in this work.
In our solution of the nonlinear FPE we can reproduce neither the available data at SPS and RHIC energies for any value of q ∈ (1, 1.5), nor the corresponding solutions from Refs. [7][8][9]. The use of three different numerical methods with coinciding outcome and various cross checks with exact analytical results ensure the accuracy of our numerical calculations.
This result casts doubt on the validity of the nonextensivity concept in statistical physics, which has often been applied to interpret observables in relativistic heavy-ion collisions. Nevertheless, formulae derived from nonextensive statistics may still be used in phenomenological fits of transverse momentum distributions in relativistic collisions because they allow to account for the observed transition from exponential to power law distributions that the e x q -function properly describes with q ∈ (1, 1.5). This transition had already earlier been modeled by Hagedorn in Ref. [26] using an equivalent, but QCD-inspired formula. The data for rapidity distributions in stopping and particle production can be described using an unmodified Fokker-Planck equation with a linear diffusion term as in Boltzmann statistics. Here we determine empirical values of the diffusion coefficient that are necessary to reproduce the measured data, thus accounting phenomenologically for both, nonequilibrium thermal processes and collective expansion, without a nonlinear diffusion term.