Nonextensive hydrodynamics of boost-invariant plasmas

We use quasiparticle anisotropic hydrodynamics to study the non-conformal and non-extensive dynamics of a system undergoing boost-invariant Bjorken expansion. To introduce nonextensivity, we use an underlying Tsallis distribution with a time-dependent nonextensivity parameter $q$. By taking moments of the quasiparticle Boltzmann equation in the relaxation-time approximation, we obtain dynamical equations which allow us to determine the time evolution of all microscopic parameters including $q$. We compare numerical solutions for bulk observables obtained using the nonextensive evolution with results obtained using quasiparticle anisotropic hydrodynamics with a Boltzmann distribution function ($q \rightarrow 1$). We show that the evolution of the temperature, pressure ratio, and scaled energy density, are quite insensitive to which distribution function is assumed. However, we find significant differences in the early-time evolution of the bulk pressure which are observed for even small deviations from the Boltzmann distribution function. Finally, we discuss the existence of non-conformal hydrodynamic attractors for the longitudinal and transverse pressures, the bulk and shear viscous corrections, and the nonextensivity parameter $q$.


I. INTRODUCTION
Ultrarelativistic heavy-ion collision experiments allow physicists to study the behavior of nuclear matter at extremely high temperatures. In such experiments, high-energy collisions of heavy nuclei are used to heat a volume of matter up to temperatures that exceed the critical temperature (T c ∼ 155 MeV) necessary to create a quark-gluon plasma (QGP) [1][2][3][4]. The bulk evolution of the QGP, as a strongly interacting state of matter, is well described by relativistic hydrodynamics [5][6][7][8][9][10][11]. In the last two decades, different frameworks have been developed to describe many heavy-ion collision observables [12][13][14][15][16][17][18][19][20] (see Ref. [21] for a recent review). One of the frameworks used is anisotropic hydrodynamics, which was introduced to take into account the fact that the QGP is a highly momentum anisotropic plasma at early times after the nuclear collision [22,23]. Motivation and the basics of anisotropic hydrodynamics are presented in detail in Ref. [24].
In recent years, 3+1D quasiparticle anisotropic hydrodynamics has been developed where three momentum-space anisotropy parameters in the underlying distribution function and a single-finite thermal mass which is fit to lattice QCD data for the equation of state are included. The results from this model have been compared with a variety of heavy-ion observables at different collision energies: Au-Au collisions at 200 GeV and Pb-Pb collisions at 2.76 TeV and 5.02 TeV [25][26][27][28][29]. In all these cases, the 3+1d quasiparticle anisotropic hydrodynamics model was able to describe the data reasonably well for many heavy-ion observables such as the spectra, mean transverse momentum of identified hadrons, multiplicities, the elliptic flow, and the HBT radii. See Refs. [30,31] for recent reviews of 3+1D quasiparticle anisotropic hydrodynamics.
One of the challenges of hydrodynamic models is to describe the spectra at intermediate transverse momentum, p T ∼ 3 GeV. As an example, in the 3+1D aHydroQP model, especially in peripheral collisions, one observes that the agreement between the aHydroQP model and experimental results is good only for p T < ∼ 1.5 GeV [29]. While the data exhibits a characteristic power-law tail, aHydroQP predictions exhibit an exponential behavior which falls faster than the experimental data does. This difference is primarily due to the fact that thermal distribution functions, which are asymptotically of exponential form, are assumed at freeze-out in aHydroQP. In contrast, when used at freeze-out, statistics based on a Tsallis distribution function can provide exponential behavior at low p T and power-law behavior at high p T [32][33][34][35].
Tsallis statistics has been used to fit the spectra more accurately over a wide range of p T in many high-energy experiments involving different systems and different collision energies, [36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51]. Rather than performing a direct fit to the data using a Tsallis distribution for the freeze-out hypersurface distribution, one would like to account for the nonextensivity in the underlying dynamical model, using Tsallis non-extensive statistics. It would, therefore, be of interest to employ the Tsallis distribution function in the dynamical model itself to describe the characteristics of the medium in both the QGP and the freeze-out phases in order to study the effects on the predicted bulk observables.
In a prior paper [52], quasiparticle anisotropic hydrodynamics was studied in 0+1D systems undergoing boost-invariant Bjorken expansion with the assumption of a Boltzmann distribution function. In this work, we derive the dynamical equations for the same systems using a Tsallis distribution function. After solving the dynamical equations, we compare the bulk observables predicted by this approach with the results presented before in Ref. [52].
We show that the evolution of temperature, pressure anisotropy, and energy density are not very sensitive to which distribution function is assumed in the aHydroQP framework. However, we find significant differences in the evolution of the bulk pressure which are observed for even slight deviations from the Boltzmann distribution function. We finally investigate the existence of dynamical attractors in both approaches. We present the evolution of the longitudinal and transverse pressures and the bulk and shear viscous stresses as a function of time in units of the local relaxation time. We find that far-from-equilibrium non-conformal attractors exist for both extensive and non-extensive statistics.
The structure of the paper is as follows. In Sec. II, we introduce the Tsallis and Boltzmann distribution functions. In Sec. III, we review the basics of the 0+1D quasiparticle anisotropic hydrodynamics model. In Sec. IV, results are presented for comparisons between the Boltzmann and Tsallis distribution functions using the 0+1D quasiparticle anisotropic hydrodynamics framework. Sec. V contains our conclusions and an outlook for the future.

II. TSALLIS DISTRIBUTION FUNCTION
The Tsallis non-extensive distribution function is given by [32,33,53] where q is the Tsallis parameter which characterizes the nonextensivity of the system. By expanding Eq. (1) around q = 1, one obtains where the leading term is the Boltzmann distribution function. Since the q → 1 limit yields the Boltzmann distribution, it is helpful to characterize the Tsallis distribution function in terms of deviations from q = 1 (which would correspond to the Boltzmann limit). For this purpose we introduce the notation δq ≡ q − 1 which is a measure of the degree of non-extensivity, so that the δq = 0 limit corresponds to the Boltzmann distribution.
Hence the Tsallis distribution can be rewritten as Henceforth all results obtained using the Tsallis distribution function shall be characterized in terms of the parameter δq, and all results must be interpreted in light of how far or near the system is from the thermal distribution. An interpretation of the Tsallis distribution and a comparison to Boltzmann distribution with illustrations of the number density are provided in Refs. [39,54,55]. For more information concerning the application of Tsallis statistics to high-energy proton-proton, proton-nucleus, and nucleus-nucleus collisions, we refer the reader to the recent article of Kapusta [56] and references therein.

III. ANISOTROPIC HYDRODYNAMICS
In the canonical anisotropic hydrodynamics approach, the distribution function is assumed to be anisotropic in momentum space i.e, in the local rest frame one has [22,23,52] f (x, p) = f eq 1 λ where λ corresponds to a temperature-like parameter, α i corresponds to the momentum anisotropy parameter in the i-th direction, and f eq is an exponential Boltzmann distribution.
In the limit where α i = 1 and λ = T , one recovers the isotropic thermal distribution. In this work, the distribution function has the same argument as in Eq. (4), but with Tsallis form rather than Boltzmann form i.e., Note that in the limit of δq → 0, Eq. (5) reproduces Eq. (4).

A. Bulk Variables
The bulk variables: number density n, energy density E, and the pressure P can be computed once the distribution function is specified. Below we list their definitions, respectively: Here dP is the Lorentz invariant momentum-space integration measure given by where N dof is the number of degrees of freedom. For compactness,Ñ ≡ N dof /(2π) 3 is used below.

B. Quasiparticle anisotropic hydrodynamics
The quasiparticle Boltzmann equation is given by [52,57,58] p where p µ is the four momentum, ∂ µ is the four derivative, (µ for spatiotemporal indices and i for spatial indices), m is the mass, f is the phase space density and the right hand side is the collision kernel containing all interactions. In this work, the collisional kernel is taken in the relaxation-time approximation (RTA), C[f ] = p µ u µ (f − f eq )/τ eq with τ eq being the relaxation time [52].
Taking the zeroth, first, and second moments of the Boltzmann equation yields, respectively, one obtains where J µ is the particle four-current, T µν is the energy-momentum tensor, and I µνλ is a rank-three tensor. They are given by In Eq. (15), a background contribution B is introduced to ensure thermodynamic consistency as explained in detail in Ref. [52]. As a result, there exists a partial differential equation relating B and the thermal mass The thermal mass m(T ) is obtained by tuning to the equation of state from lattice QCD calculations [59]. The procedure for performing this extraction is explained in detail in Ref. [52].

C. 0+1D Quasiparticle anisotropic hydrodynamics
In this work, we will limit ourselves to boost-invariant systems (0+1D). In this case, the energy density, transverse pressure, and longitudinal pressure are given by [52] E =H 3 (α,m) λ 4 + B , In 0+1D, Eq. (17) can be written as For definitions of the various H-functions appearing above, we refer the reader to App. A.
In the case where a Tsallis distribution is assumed, the q dependence is implicit in f (x, q), which the H-functions are based upon.

D. The dynamical equations in 0+1D
In the 0+1D case, there are five dynamical variables λ, T , α x , α z , and q. Hence we need 5 dynamical equations which can be obtained from the zeroth, first and second moments of the Boltzmann equation.
From the zeroth moment, one obtains, where n eq = 4πÑ T 3m2 eq K 2 (m eq ) . and From the first moment, the conservation of energy-momentum ∂ µ T µν = 0, one obtains, From the second moment, Eq. (13), one obtains where with K 3 is a modified Bessel function of the second kind and I is given by When f (x) is in the Tsallis form, the integral in Eq. (30) can be done numerically, where, upon using the Boltzmann distribution function, the following result is obtained The fifth equation necessary can be obtained by using the matching condition which reflects energy-momentum conservation: By expanding Eqs. (20), (24), (25), (26) and (32), the final equations can be written in the following format: where τ eq is a function of T (τ ) and the temperature-dependent quasiparticle mass m(T ) [52].
When using a Boltzmann distribution, there are four dynamical equations obtained solely from the first and second moments; derived in detail in Ref. [52]. They also can be obtained from Eqs. (34)-(37) by setting ∂ τ q = 0 and using the Boltzmann distribution in the definitions of theH, Ω, and I functions. The reader may refer to App. A for more details about the special functions specified above.
One should note that, from the definitions of the bulk variables, one needs to obtain B in order to obtain the full energy density and pressures. For this purpose, we integrate the dynamical equations to very late times (τ f = 100 fm/c) and then integrate Eq. (19) backwards in time using B(τ f ) = B eq (T (τ f )) since B eq (T = 0) = 0 [52].

IV. RESULTS
In this section, we present comparisons of the 0+1D quasiparticle anisotropic hydrodynamics model using both the Boltzmann and Tsallis distribution functions. For reference, we compare the results of this work to the results presented in Ref. [52], where a Boltzmann distribution function is used instead of the Tsallis distribution function. We first solve the dynamical equations presented in the last section subject to this set of initial conditions: T 0 (τ 0 ) = 600 MeV, α x (τ 0 ) = 1, and α z (τ 0 ) = 1 at an initial proper time τ 0 = 0.25 fm/c.
Then we show differences in the temporal evolution of temperature, energy density, pressure anisotropy, and bulk pressure, which is defined as pressure, and the transverse anisotropy parameter α x . In this figure, the shear viscosity to entropy density ratio is taken to be 4πη/s = 1. As can be seen from the figure, the effective temperature evolution in both approaches is identical whereas some differences are seen in the evolution of the scaled energy density and pressure anisotropy. The maximum difference for the scaled energy density is roughly ∼ 1% at late times whereas the maximum difference for the pressure anisotropy is roughly ∼ 10% at τ ∼ 0.5 fm/c. The temporal evolution of the bulk pressure shows clear differences between the two methods which means the bulk pressure is sensitive to which distribution function is assumed. We note that differences in the bulk evolution exist between these two approaches even for very small initial q values such as δq = 0.001 (not shown here). As shown in Ref. [60], differences in the evolution of the bulk pressure can have a direct impact on the primordial particle spectra and hence this may be important phenomenologically. Next, we show that, at late times, both the Tsallis parameter q and the transverse anisotropy parameter α x approach unity as expected by using the RTA approximation.
In Fig. 2, we show the proper-time evolution of the same quantities as in Fig. 1 with the shear viscosity to entropy density ratio is taken to be 4πη/s = 5. As can be seen from this figure, deviations between the two approaches become more pronounced when increasing the shear viscosity to entropy density ratio. We note, as can be seen from panels (e) and (f), that it takes a longer time for both approaches to converge to a universal results, compared to results in Fig. 1 for the same panels.
We note that in the Tsallis approach, there are numerical instabilities appearing at late times resulting from integrals that depend on derivatives of f T as shown in App. A. Although they do not affect the late time behavior of some quantities such as the pressure anisotropy, they affect late time behavior of the bulk pressure. To control this numerical issue, we take two steps. First, we expand f T (x, q) around q → 1 up to the sixth-order and use this approximated f T when q is close to 1, this point is taken to be q c = 1.00001. Second, since the integrands are practically zero at large x, we put a limit on the integrations instead of integrating to ∞. With these changes, we still see some numerical issues at late times, however, we always terminate evolution before they become too serious.

B. Non-equilibrium attractors
We next turn to a discussion of non-equilibrium attractors. For this purpose we will make plots of the shear viscous stress, which is defined as In Figs. 3-4, we investigate the attractor behavior corresponding to a variety of initial conditions using Boltzmann and Tsallis distributions, respectively. In both figures, we assumed 4πη/s = 3 and δq(τ 0 ) = 0.10. We show the scaled time evolution τ /τ eq and for reference, in some cases, we also show τ T evolution which is the limit for the conformal relaxation time, i.e. τ eq ∼ 1/T (τ ). We note that the black-solid line is initially isotropic, i.e. α x = α z = 1. As a result, P L /P eq = P T /P eq = 1, whereas the bulk and shear viscous stress are both equal to zero. In the same figures, the other colored lines correspond to different initial conditions where α x = α z = 1. Finally, to put the numbers presented in perspective, we note that at τ 0 , τ T = 0.15 and τ /τ eq = 0.628217 and, at very late times, e.g. τ = 100 fm/c, τ T ∼ 15 and τ /τ eq ∼ 35.
In Fig. 3, top and middle rows, we show the evolution of the scaled longitudinal pressure P L /P eq and the scaled transverse pressure P T /P eq , respectively. Both P L /P eq and P T /P eq Top row: Evolution of P L /P eq as a function of τ /τ eq and τ T , in the left and right panels, respectively. Middle row: Evolution of P T /P eq as a function of τ /τ eq and τ T , in the left and right panels, respectively. Bottom row: Evolution of Π/P eq and π/P eq as a function of τ /τ eq , in the left and right panels, respectively.
are shown as a function of τ /τ eq and τ T in the left and right panels, respectively. We see that all results for different set of initial conditions converge to a universal curve at early times τ ∼ 4 fm/c. In the bottom row, we plot results of the scaled bulk pressure Π/P eq Evolution of q and P T /P as a function of scaled timed τ /τ eq , left and right panels, respectively. Note that q 0 is different in each curve as can be seen from left panel whereas other initial conditions are the same as the ones used in Fig. 4. and the shear stress π/P eq as a function of τ /τ eq . We note that results for π/P eq approach a universal attractor at an earlier rescaled time τ /τ eq ∼ 4 while the Π/P eq converges later, only for τ /τ eq > 10 due to the system reaching isotropic equilibrium. In the context of initial conditions used here τ /τ eq = 4 corresponds to τ ∼ 4 fm/c, whereas τ /τ eq = 10 corresponds to τ ∼ 20 fm/c.
In Fig. 4, we use a Tsallis distribution and present the same plots as shown in Fig. 3.
From all panels, we see the existence of a universal attractor even for this nonextensive farfrom equilibrium approach. The late time differences in the Π/P eq evolution seen in panel (e) could be purely numerical in origin; however, the differences seen at earlier rescaled times indicate that, strictly speaking, there does not exist a hydrodynamic attractor for the evolution of the bulk viscous pressure. Despite this, similar to as was found using exact solutions in the relaxation time approximation [61], we find that the transverse and longitudinal pressure ratios shown in panels (a) and (c) and the shear correction shown in panel (f) suggest the existence of a non-equilibrium hydrodynamic attractor even when the system has a realistic non-conformal equation of state. This is true for both Boltzmann and Tsallis statistics.
Finally, in Fig. 5, we change the initial condition for the Tsallis parameter, q 0 , for each curve, with the other initial conditions being the same as in Fig. 4. In the left panel, we show the scaled time evolution of the Tsallis parameter q. As can be seen from this figure, all results converges quickly to unity at τ /τ eq ∼ 5. In the right panel, as a cross check, we plot the scaled time evolution of P T /P eq for this set of initial condition and we can see that a universal attractor still exists. Figure 5 demonstrates that an attractor exists for the non-extensivity parameter q, with q approaching unity from above. This is a non-trivial finding and suggests that one can constrain the late-time value of q used in phenomenological applications without a full dynamical simulation.
We note that fits to the experimental data indicate that the fitted Tsallis parameter q is always close to unity meaning δq is small. For example, in Au+Au collisions at 200 GeV, fits to the spectra result in δq = 0.015 and δq = 0.086 at 10-20% and 60-80%, respectively [40]. In Pb-Pb collisions at 2.76 TeV, fits to the spectra result in δq = 0.1363 at 5-10 % centrality [62].
Moreover, in p-p collisions at 5.02 TeV and 13 TeV, fits to the spectra data result in δq ∼ 0.12 which increases as a function of multiplicity to reach δq ∼ 0.14 [47]. For results presented in this work, we used different values of δq ∈ {0.15, 0.10, 0.08, 0.06, 0.04, 0.02}. In practice, one would like to go to higher values of δq; however, for large δq the moment expansion of the distribution is ill-defined due to UV divergences. In addition, we find that our numerics become quite unstable for large δq close to this limit. This is expected since the H functions introduced in the App. A, which are part of the dynamical equations, each have their own convergence intervals. The convergence intervals for the bulk variables in rather simplified systems are shown in Ref. [53]. As an example, in the massless limit, the number density is convergent only in the interval 0 ≤ q ≤ 1.5.

V. CONCLUSIONS AND OUTLOOK
In this paper, we reviewed the basics of quasiparticle anisotropic hydrodynamics which is based on the self-consistent introduction of a single temperature-dependent quasiparticle mass for the degrees of freedom in the Boltzmann equation. The temperature dependence of the quasiparticle mass was determined uniquely by matching to lattice QCD calculations of the equation of state. We then used 0+1D quasiparticle anisotropic hydrodynamics as a dynamical model, with the underlying distribution given by either a Tsallis distribution function or its δq → 0 limit, which corresponds to a Boltzmann distribution function. We then compared the temporal evolution of temperature, pressure anisotropy, energy density, and bulk pressure predicted by the two approaches at different shear viscosity to entropy density ratios. This work demonstrates that the temperature, pressure anisotropy, and energy density evolutions are not very sensitive to which distribution function is used in the model. However, we found that the bulk evolution is sensitive to which distribution function is assumed even for very small δq ∼ 0.001.
In the last section of the paper we demonstrated the existence of hydrodynamic attractors in both non-conformal cases. We found that early-time hydrodynamic attractors exist for the scaled longitudinal and transverse pressures even though the system is non-conformal and non-extensive. We also found that, although the bulk viscous pressure does not have an early-time hydrodynamic attractor, the shear stress seems to converge to an early-time hydrodynamic attractor. These results are in agreement with earlier findings which made use of exact solutions to the 0+1D RTA Boltzmann equation with Boltzmann statistics [61].
Finally, by varying q 0 , we found that there exists a hydrodynamic attractor for the Tsallis parameter q when q is plotted as a function of rescaled time. This observation can help to strongly constrain late-time values of q for phenomenological applications.
As shown in Refs. [52,60], differences in the evolution of the bulk pressure can have a direct impact on the primordial hadron spectra. Looking forward, these differences could allow a determination of the optimal form of the distribution function on the freeze-out hypersurface. This could be done by modifying the existing 3+1d aHydroQP code where the distribution function is assumed to be in the Boltzmann form [18], to study the effect of Tsallis statistics on heavy-ion observables such as the spectra especially at intermediate p T ∼ 3 GeV. On the theory front, it would be quite interesting to study the conformal hydrodynamic attractor in the case of Tsallis distribution at both zero and finite chemical potential as performed in Refs. [63][64][65]. Additionally, one also may look for the existence of the exact solutions to the Boltzmann equation using a Tsallis distribution function at small δq [66]. These projects are planned for future follow-up work. Finally, we note that in deriving the dynamical equations, the following identities are needed: ∂H 2 (y, z) ∂y = 1 y H 2 (y, z) + H 2L (y, z) , (A18) and

Second moment
The special functions related to the second moment of the Boltzmann equation are I (m) ≡m ∂mI =m 2 dpp 4 p 2 +m 2 f p 2 +m 2 , and I q (m) ≡ ∂ q I = dpp 4 f q p 2 +m 2 . (A24)