New Exact Solutions of Relativistic Hydrodynamics for Longitudinally Expanding Fireballs

Abstract

We present new, exact, finite solutions of relativistic hydrodynamics for longitudinally expanding fireballs for arbitrary constant value of the speed of sound. These new solutions generalize earlier, longitudinally finite, exact solutions, from an unrealistic to a reasonable equation of state, characterized by a temperature independent (average) value of the speed of sound. Observables such as the rapidity density and the pseudorapidity density are evaluated analytically, resulting in simple and easy to fit formulae that can be matched to the high energy proton–proton and heavy ion collision data at RHIC and LHC. In the longitudinally boost-invariant limit, these new solutions approach the Hwa–Bjorken solution and the corresponding rapidity distributions approach a rapidity plateaux.

1. Introduction

Some of the most renowned theoretical papers in high energy heavy ion physics deal with exact solutions of perfect fluid hydrodynamics for a 1 + 1 dimensional, longitudinally expanding fireball. In high energy collisions involving strong interactions, statistical particle production rates were noted by Fermi already in 1950 [1]. Landau, Khalatnikov and Belenkij predicted as early as in 1953–1956 [2,3,4] that, in these collisions, not only global but also local thermal equilibrium will be a relevant concept and the related perfect fluid hydrodynamical modelling will provide the framework for the future analysis of experimental data.

Landau’s prediction about the perfect fluid behaviour in high energy heavy ion collisions was, at that time quite unexpectedly, fully confirmed by the discoveries after the first four years of data-taking at Brookhaven National Laboratory’s Relativistic Heavy Ion Collider (RHIC) where the picture of a nearly perfect strongly interacting fluid emerged from the first years of observations of the four RHIC collaborations, BRAHMS [5], PHENIX [6], PHOBOS [7], and STAR [8].

When the experiment called “A Large Ion Collider Experiment” (ALICE) Collaboration reported first elliptic flow measurement at the Large Hadron Collider (LHC) in Pb + Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV [9], they also noted the similarity of the transverse momentum dependence of the elliptic flow to earlier, lower energy data at RHIC and noted its consistency with predictions of hydrodynamic models, confirming the creation of a nearly perfect fluid of strongly interacting Quark Gluon Plasma (sQGP) also at higher energies. More detailed measurements by ALICE, and the experiments called “A Toroidal LHC ApparatuS” or ATLAS and “Compact Muon Solenoid” (CMS) [10,11,12] gave further support to the hydrodynamical paradigm in high energy heavy ion collisions at RHIC and LHC. The domain of the validity of this paradigm is currently being extended to describe proton–proton, as well as to proton/d/${}^3$He-nucleus collisions, from the highest LHC energies [12] down to the lowest RHIC energies [13,14,15].

The status of the applications of hydrodynamical modelling to high energy collisions was reviewed recently [16], also including a review of exact solutions, but focussing mostly on the numerical solutions of relativistic hydrodynamics, as well as on the conceptual and open questions.

To define the context of our new exact solutions described in the body of the manuscript, let us dedicate the next section to a brief, but specialized overview of the status of the search for exact solutions of fireball hydrodynamics and the application of these results to high energy hadron-proton and heavy ion collisions. The equations of perfect fluid hydrodynamics are summarized in the subsequent section. Next, we present our new exact, accelerating family of solutions of relativistic hydrodynamics, pointing out their finite range in rapidity, but without detailing their derivation. Then, we evaluate the observable rapidity and pseudorapidity densities near midrapidity, and discuss and compare the results with experimental data on proton–proton collisions at the $\sqrt{s}=$7 and 8 TeV at LHC. Finally we summarize and conclude.

2. A Brief Overview of the Status of the Field

2.1. Exact Solutions with Boost-Invariant Flow Profiles

Hydrodynamics typically predicts strong anisotropies in the angular distribution of produced particles, with a predominantly longitudinal expansion, giving a motivation to study exact solutions of perfect fluid hydrodynamics in 1 + 1 dimensions, characterized by a temporal and a spatial dimension, where the spatial dimension corresponds to the main axis of the expansion, usually taken as the longitudinal direction.

The ultra-relativistic limit of the hydrodynamical models was, as far as we know, solved for the first time by Hwa as early as in 1974. Hwa noted that, in the limit of infinite colliding energies, the resulting rapidity distribution of produced hadrons will be flat, a constant value, also called the rapidity plateaux [17]. Independently, but several years later, Bjorken rediscovered the same boost-invariant asymptotic solution [18]. In addition, rather importantly, Bjorken also realized that the rapidity density can be used to estimate the initial energy density in these reactions [18]. By now, the resulting Bjorken energy density estimate became the top-cited phenomenological formula in high energy heavy ion physics.

The Hwa–Bjorken boost-invariant solution was extended to include transverse expansion dynamics as well. Certain exact solutions were obtained by Bíró that considered the case of transverse expansion on top of a longitudinally boost invariant flow during a time period with constant pressure. Such a scenario was thought to be relevant during the course of a first order QGP–hadron gas phase transition [19]. Gubser found a similar, longitudinally boost-invariant but transversally expanding solution for a massless relativistic gas, characterized by the equation of state $\epsilon =3p$, where $\epsilon$ denotes the energy density and p stands for the pressure. This choice of the equation of state corresponds to conformally symmetric solutions [20,21], where the transverse flow is initially vanishing, and it is increasing at smaller transverse distances nearly linearly, that however it does not keep on increasing but changes to decreasing to zero at large distances. Solutions for very viscous plasmas (cold plasma limit) were also found that recover the perfect fluid solution of Gubser in the limit of vanishing viscosity [22].

Using the equation of state of a relativistic ideal gas and conformal symmetry, analytic axially symmetric solutions with non-zero vorticity was described in Ref. [23] and detailed in Ref. [24], with the aim that the phenomenological extension of these solutions to longitudinally boost-invariant but elliptically expanding case may give clues to understand vorticity effects and the elliptic flow data in heavy ion collisions at RHIC and at LHC. Results from this approach were reviewed briefly in Ref. [25].

A new analytic solution of relativistic hydrodynamics for perfect fluids with conformal symmetry and equation of state for a massless relativistic gas, $\epsilon =3p$ was presented very recently [26]. This new solution corresponds to rotations of the fluid around two different vortices (axis), generalizing simultaneously the non-rotating solutions of Refs. [20,21], as well as the rotating solutions described in [27] and more recently described and detailed in Ref. [23,24]. The possibility of describing relativistic fluids with one and two vortices prompts the question: Can generalized exact formulae can be derived that describe expanding and locally rotating relativistic fluids that may have more than two vortices?

In this context, it is also of importance to mention that some of the exact solutions include solutions with the lattice QCD equation of state, using actually an even broader class of equations of state, where the speed of sound as well as the constant of proportionality between the pressure and the energy density may be an arbitrary but temperature dependent function. As far as we know, the first solution with a temperature dependent speed of sound was given in [28], using the equation of state $\epsilon =\kappa \left(T\right)p$ for a system that had conserved charge density n and the pressure was given as $p=nT$. Independently, this kind of equation of state was generalized for the ${\mu }_B=0$ case, when only the entropy density $\sigma$ is locally conserved, and the pressure is proportional to the entropy density, $p\propto \sigma T$ with $\epsilon =\kappa \left(T\right)p$ in ref. [29], opening the way to use lattice QCD equation of state in several subsequent exact solutions of relativistic as well as non-relativistic hydrodynamics [30,31,32].

Actually, the non-vanishing total angular momentum of the fireball created in high energy heavy ion collisions is an interesting problem in itself. It is clear that infinite, longitudinally boost-invariant sources cannot have an angular momentum perpendicular to the impact parameter plane, because their moment of inertia in the beam direction is infinite. Thus, the search for rotating solutions is indirectly coupled to the search for non-boost-invariant, finite exact solutions of fireball hydrodynamics.

The first, rotating exact solution of relativistic hydrodynamics was found by Nagy, when searching for simultaneous solutions of relativistic hydrodynamics and the relativistic collisionless Boltzmann equation [27].This searching method generalized the non-relativistic search method of Csizmadia, Csörgő and Lukács [33] to the relativistic kinematic domain, with surprising success.

2.2. Recent Results on Rotating, Non-Relativistic Solutions

Several exact solutions of rotating, non-relativistic finite fireballs were obtained recently, both for spheroidal and for ellipsoidal symmetries, corresponding to expanding and rotating spheroidal [30] or triaxial, ellipsoidal density profiles [31], which all go back to the spherically symmetric, non-rotating and collisionless solution of fireball hydrodynamics of Ref. [33]. These solutions are finite and aspire to explain some of the scaling properties of the single particle slope parameters [32,33], the elliptic flow and the Hanbury Brown–Twiss (HBT) radii.

Recently, new exact solutions including a realistic, lattice QCD based equation of state and the effects of angular momentum [34] were found to describe well the mass systematics of single particle spectra from a strongly interacting quark-gluon plasma being converted to a multi-component hadron gas [32]. It is interesting to note that, although these non-relativistic exact solutions of rotating fireball hydrodynamics were found very recently, they were obtained using analogies with relativistic exact solutions of rotating fireball hydrodynamics, that were found first by Nagy in 2009, when looking for simultaneous solutions of relativistic hydrodynamics and the relativistic collisionless Boltzmann equation [27].

As the medium cools down, quark and gluon degrees of freedom get confined to hadrons. After such a re-confinement, soon a hadrochemical freeze-out may be reached, where the equations describing local energy–momentum and entropy conservation are supplemented with the continuity equation for each of the locally conserved particle densities. This problem was considered recently in Refs. [32,34], both in the relativistic and in the non-relativistic kinematic case, successfully explaining the mass-dependence of the slope parameters of single-particle spectra from an exact solution of non-relativistic hydrodynamics, showing also that many features of these solutions are similar to the system of equations of the relativistic fluid dynamics.

One of the important observations is that, at late time, these non-relativistic solutions approach an asymptotic, three-dimensional Hubble or rotating Hubble flow profiles, which are solutions of hydrodynamics in the relativistic kinematic domain as well. In the longitudinal direction, profiles of the Hubble flow and that of the Hwa-Bjorken solution coincide. One of the important limitations of the boost-invariant Hwa–Bjorken solution [17,18] is its infinite range along the axis of expansion. In contrast, the non-relativistic, spherically symmetric Zimányi–Bondorf–Garpman solution [35] is similarly simple, but it is a compact and finite, spherically symmetric solution. Actually, the velocity field is a linear, Hubble-type flow field with $\mathbf{v}\propto \mathbf{r}$ in both of the relativistic Hwa–Bjorken and in the non-relativistic Zimányi–Bondorf–Garpman solution [17,18,35]. Thus, it would be worthwhile to search for the relativistic generalizations of the non-relativistic but rotating exact solutions given in Refs. [30,31,32,34].

2.3. Relativistic Solutions without Longitudinal Boost-Invariance and Rotation

One of the natural themes when searching for new exact solutions of relativistic hydrodynamics was to generalize the Hwa–Bjorken boost-invariant longitudinal flow to a finite, realistic, accelerating longitudinal flow profile. Initially, such solutions were found by breaking the boost-invariance of the temperature or the (entropy)density distributions, while keeping the pressure and the flow-field boost-invariant [36,37,38].

In these solutions, the Hubble flow field remained an important element, but one obtained generalized, spatially inhomogeneous temperature and density profiles in new, self-similar, 1$+$1 dimensional and 1$+$3 dimensional axially symmetric solutions by Csörgő, Grassi, Hama and Kodama in Refs. [36,37], or to 1$+$3 dimensional ellipsoidally symmetric solutions by Csörgő, Csernai, Hama and Kodama of Ref. [38]. Below, we refer to these as the CGHK and the CCHK solutions, respectively. In these solutions, the equation of state was reasonable as well, as the constant of proportionality between the pressure and the energy density had arbitrary value. The transverse momentum spectra of identified pions, kaons and protons, as well as elliptic flow data and HBT radii, the characteristic length-scales measured by Bose–Einstein correlation functions were shown to be well described by the CCHK solution [38], as demonstrated in [39]. Even the direct photon spectra and the elliptic flow of direct photons was shown to be well described by this CCHK solution [40]. However, from the theoretical point of view, the main limitation of the CCHK solution was its non-accelerating property, due to the presence the Hubble flow field, that can be seen as an asymptotic attractor for late stages of fireball explosions.

Let us note here that hydrodynamics was successfully applied to describe the double-differential rapidity and transverse mass spectra in hadron + proton collisions at the surprisingly low energy of $\sqrt{s}\simeq 22$ GeV already in 1997, using the Buda–Lund hydro model [41], which also predicted a special coupling between longitudinal dynamics and transverse expansion. The hydrodynamically predicted decrease of the slope of the transverse momentum spectra at forward and backward rapidities was experimentally observed [42]. Although phenomenologically very successful, this Buda–Lund hydrodynamical model used a boost-invariant longitudinal flow profile, but a non-boost-invariant longitudinal chemical potential or density. Subsequently, the time evolution, corresponding to the Buda–Lund hydro model of Ref. [41], was related to the time evolution in the CCHK exact solution [38].

2.3.1. Landau Hydrodynamics

One of the research directions to describe the finite rapidity or pseudorapidity distributions was due to Wong and collaborators, who first improved the evaluation of the rapidity density [43] and subsequently studied the matching of the inwards-moving shock-wave of Khalatnikov’s solution with the outward-moving regular solution used by Landau to obtain an approximately Gaussian rapidity distribution, using the massless relativistic ideal gas equation of state, $\epsilon =3p$ [44].

By using the revised Landau hydrodynamic model and taking into account the effect of leading particles as independent Gaussian sources at forward and backward rapidities, Jiang, et al. [45] explained the pseudorapidity distributions of produced charged particles at different centralities in $\sqrt{s_{NN}}=200$ GeV Cu + Cu and Au + Au collisions at RHIC. Similarly good descriptions of the pseudorapidity distributions were obtained at lower colliding energies, $\sqrt{s_{NN}}=130$ and 62 GeV [46].

Starting from about 2014, the Landau hydrodynamics was revisited and tested successfully on proton–proton collisions as well, from $\sqrt{s}=$ 23 to 900 GeV [47,48]. The low energy limit of these good fits also moved further down in the case of heavy ion collisions to $\sqrt{s_{NN}}=$ 19.6 and 22.4 GeV Au + Au and Cu + Cu collisions [49], and the range also was extended upwards to $\sqrt{s_{NN}}=$ 2.76 TeV Pb + Pb collisions at CERN LHC. To achieve a successful fit, a Gaussian contribution explained as leading particle contribution had to be added to the Landau hydrodynamical calculations both at forward and at backwards rapidities [50]. Landau hydrodynamics was generalized also for an arbitrary constant speed of sound and it was shown to give a similarly reasonable description of the pseudorapidity densities at 2.76 TeV Pb + Pb collisions [51].

2.3.2. Results from the BJP Solution

Another accelerating, one-parameter family of analytic solutions interpolating between the boost-invariant Bjorken picture and the non boost-invariant one, similar to the Landau profile was described by Bialas, Janik and Peschanski in [52], referred to as the BJP solution. Using the Khalatnikov potential, an analytic formula was derived from the BJP solution for the rapidity dependence of the entropy density, $dS/dy$ [53].

Recently, it was found that the BJP solution describes the longitudinal evolution in relativistic heavy ion collisions both at RHIC and at LHC, when a leading particle effect is included, for a realistic, constant value of the speed of sound, including the centrality dependence in Cu + Cu and Au + Au collisions at $\sqrt{s_{NN}}=200$ GeV and in Pb + Pb collisions at $\sqrt{s_{NN}}=2.76$ TeV [54]. The same BJP family of solutions were shown to describe the longitudinal evolution and the resulting pseudorapidity distributions in proton–proton collisions as well, from 23 GeV to 7 TeV [55], when the leading particle effect was corrected for in the forwards and backwards direction. Thus, the BJP solution of relativistic hydrodynamics describes well the longitudinal observables both in proton–proton and in heavy ion collisions, from the lowest to the highest presently available energies [56,57].

The elliptic and higher order flows were also obtained analytically from the BJP family of solutions, using the assumption that the entropy is transversally conserved. This assumption holds for the asymptotic Hwa–Bjorken flows but, as far as we can see, it is not a generally valid property of accelerating longitudinal flows, as it neglects the important correction due to work done by central fluid elements on the surface [58]. Despite these successes of the BJP solution, these solutions have not yet been connected to possible estimations of the initial energy density of high energy proton–proton or heavy ion collisions.

2.3.3. Results from the CNC Solution

Motivated by searching for solutions with finite rapidity distributions and corrections to the initial energy density estimate of Bjorken, an exact and explicit, longitudinally accelerating solutions of relativistic hydrodynamics, a one-parameter family of analytic solutions was presented for 1$+$1 dimensional explosions by Csörgő, Nagy and Csanád in ref. [59] and detailed in Ref. [60], referred to as the CNC solution. This CNC solution generalized the Bjorken flow field of $v_z=tanh\left({\eta }_x\right)=r_z/t$ to $v_z=tanh\left(\lambda {\eta }_x\right)$, where $v_z$ is the longitudinal component of the four-velocity field, ${\eta }_x=(1/2)ln\left[(t+r_z)/(t-r_z)\right]$ is the space-time rapidity and $\lambda$ is a parameter that characterizes the acceleration of the fluid. This one parameter family of solutions had five different domains of applicability, where different domains were characterized by the value of the acceleration parameter $\lambda$, the number of spatial dimensions d and the parameter of the equation of state, $\kappa$.

Almost simultaneously, Borshch and Zhdanov published exact solutions using the superhard $\kappa =1$ equation of state that included the CNC solution [59,60] in certain limiting cases, and some other more general solutions also where the equation of state was more realistic [61], but the phenomenology related to observables in high energy heavy ion collisions has not yet been developed as far as we know. It would be interesting to see how these Borshch–Zhdanov (BZ) solutions compare to experimental data at RHIC and at LHC.

It is also interesting to note that the BZ solutions also were given in a form that is well suited to study fluctuations or ripples, sound waves propagating on the top of an already known solution. The 1+1 dimensional problem of relativistic hydrodynamics was solved for arbitrary intitial conditions in the $\kappa =1$ case both in the CNC [59,60] and in the BZ papers [61]. Ripples on known 1+3 dimensionally expanding solutions of relativistic hydrodynamics were studied recently in refs. [62,63,64].

In one out of these five classes of the CNC solutions [59], the parameter $\lambda$ was arbitrary; its value could be determined from fits to experimental data. This CNC solution not only resulted in a realistic and straightforwardly usable result for both the rapidity and the pseudorapidity densities, but these results were also applied to get important and large corrections (factors of 5–10) to the initial energy density estimate of Bjorken. These corrections were argued to take into account the work done by the central fluid elements on the surface. Unfortunately, these CNC solutions have also a major shortcoming, namely that the acceleration parameter $\lambda$ became a free fit parameter only for the superhard equation of state of $\kappa =1$, $\epsilon =p$. In this case, the speed of sound is equal to the speed of light, so the investigation was thought to be rather academic. Nevertheless, the equation of state dependence of the initial energy density estimate was exactly determined for $\kappa =1$ case, and it was determined to yield important corrections to Bjorken’s energy density estimate, a factor of 10 at RHIC energies. However, the dependence of this correction factor on the equation of state has not yet been determined exactly, this factor was only conjectured so far, but its numerical value was found to yield large corrections, of the order of 15 [65,66].

In this work, we present a solution that may be the basis of an exact derivation of an equation of state and acceleration dependent correction to Bjorken’s initial energy density estimate, as one of the motivations of our study was to rigorously derive the equation of state dependence of the corrections to Bjorken’s energy density estimate. Thus, the search for accelerating solutions has been continued, and some of the new results that have been achieved so far are presented in the body of this manuscript.

Recently, the formulas obtained by the CNC solution, that include the work effect but have an unrealistic, superhard equation of state, were also shown to describe surprisingly well the pseudorapidity densities in proton–proton collisions at the LHC energies [66,67]. At the same time, the pseudorapidity densities of heavy ion collisions were described well, including RHIC and LHC energies as well, from $\sqrt{s_{NN}}=130$ GeV to $\sqrt{s}=8$ TeV. These results may indicate, that the longitudinal expansion dynamics of high energy proton–proton and heavy ion collisions is surprisingly similar. In addition, they indicate that the CNC solution works surprisingly well, much better than expected when compared to the experimental data, so in some sense its final stage at the freeze-out is likely very similar to a final state obtained from solutions of relativistic hydrodynamics with more realistic equations of state. Thus, the success of the CNC fits of Ref. [67] motivated our current paper to search for accelerating, finite, exact solutions of 1$+$1 dimensional relativistic hydrodynamics with reasonable (but yet temperature independent) value for the speed of sound.

Motivated by the success of Landau hydrodynamics [44,45] and pointing out the surprising and almost unreasonable success of the CNC solution [59,60,65] to describe the pseudorapidity distributions from proton–proton to heavy ion solutions [66,67], even with a superhard equation of state, $\epsilon =p$, we suspected that the CNC solution may be very close to certain analytic solutions that have arbitrary equation of state but the fluid rapidity $\mathrm{\Omega }$ is still nearly proportional to the coordinate-space rapidity ${\eta }_x$. Indeed, we found that such exact solutions exist and they are surprisingly close in shape to the CNC solutions, given the approximate proportionality between the fluid rapidity and the space-time rapidity in these solutions, if we determine these not too far away from mid-rapidity.

3. Equations of Relativistic Hydrodynamics

For perfect fluids, relativistic hydrodynamics expresses the local conservation of energy, momentum and entropy:

$$\begin{array}{cc}\hfill {\partial }_{\nu }{T}^{\mu \nu }& =0,\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill {\partial }_{\mu }\left(\sigma u^{\mu }\right)& =0,\hfill \end{array}$$

(2)

where $T^{\mu \nu }$ is the energy–momentum four-tensor, $\sigma =\sigma \left(x\right)$ is the entropy density and $u^{\mu }$ stands for the four-velocity normalized as $u^{\mu }{u}_{\mu }=1$. These fields depend on the four-coordinate $x^{\mu }=(t,\mathbf{r})=\left(t,r_x,r_y,r_z\right)$. The four-momentum is denoted by $p^{\mu }=(E_p,\mathbf{p})=\left(E_p,p_x,p_y,p_z\right)$ with an on-shell energy $E_p=\sqrt{m^2+{\mathbf{p}}^2}$, where m is the mass of a given type of observable particle. Equation (1) stands for local conservation of energy and momentum, while Equation (2) expresses the local conservation of entropy density, and the lack of dissipative terms in perfect fluid hydrodynamics.

The energy–momentum four-tensor of perfect fluids is denoted by $T^{\mu \nu }$ which reads as follows:

$$T^{\mu \nu }=\left(\epsilon +p\right)u^{\mu }{u}^{\nu }-p{g}^{\mu \nu },$$

(3)

where the metric tensor is denoted by $g^{\mu \nu }=\mathrm{diag}(1,-1,-1,-1)$, and as already explained in Section 2, the energy density is denoted by $\epsilon$ and p stands for the pressure.

The basic equations of relativistic, perfect fluid hydrodynamics is frequently applied also in the case when there are (perhaps several) locally conserved charges. For example, after hadrochemical freeze-out at the temperature $T_{chem}$, the particle number density $n_i\left(x\right)$ is locally conserved for each of the frozen-out hadronic species indexed by i. In this case, the entropy equation is supplemented by the local charge conservation equations, expressed as follows:

$${\partial }_{\mu }\left(n_i{u}^{\mu }\right)=0,\mathrm{for}T\le T_{chem},i=1,2,\dots ,j$$

(4)

where $n_i$ is the particle density of the i-th hadron, and j counts that how many kind of hadrons are frozen out hadrochemically. For such a mixture of hadrons, $m_i$ denotes the mass for hadron type i.

The equation expressing the local conservation of energy and momentum can be projected to a component parallel to $u^{\mu }$ that yields the energy equation:

$${\partial }_{\mu }\left(\epsilon u^{\mu }\right)+p{\partial }_{\mu }{u}^{\mu }=0,$$

(5)

while the component pseudo-orthogonal to the four-velocity field yields the relativistic Euler equation:

$$\left(\epsilon +p\right)u^{\nu }{\partial }_{\nu }{u}^{\mu }=\left(g^{\mu \nu }-u^{\mu }{u}^{\nu }\right){\partial }_{\nu }p.$$

(6)

We have five differential equations for six independent quantities, the three spatial components of the velocity field and the entropy density, the energy density and the pressure, ($v_x\left(x\right)$, $v_y\left(x\right)$, $v_z\left(x\right)$, $\sigma \left(x\right)$, $\epsilon \left(x\right)$, and $p\left(x\right)$). Thus, the system of partial differential equations that corresponds to relativistic hydrodynamics is closed by the equation of state (EoS), that defines the connection between two of the six unknown fields. Throughout this paper, let us assume that the energy density is proportional to the pressure:

$$\epsilon =\kappa p,$$

(7)

where $\kappa$ is assumed to be a temperature independent constant. Thus, in this work, we search for 1$+$1 dimensional exact solutions of relativistic hydrodynamics that correspond to a generalized, $\epsilon =\kappa p$ equation of state, where the realistic value of the speed of sound is about $c_s=1/\sqrt{\kappa }=0.35\pm 0.05$ [68], so the reasonable range of $\kappa$ is approximately from 6 to 11. This EoS closes the system of partial differential equations that constitute relativistic hydrodynamics.

In a fully realistic solution, we should use the lattice QCD equation of state, with the speed of sound being a temperature dependent function, similar to Refs. [28,29,30,31,32,34]. However, in the current work we approximate $c_s\left(T\right)$ with a temperature independent constant $c_s$. We postpone the analysis of the case of $c_s\left(T\right)$ for a later, more detailed investigation. The success of the BJP [54,55,56] and the CNC solutions [66,67] in describing pseudorapidity densities using a temperature independent constant of proportionality $\kappa$ between the energy density and the pressure provides further, independent support for these investigations.

With the help of the fundamental equation of thermodynamics, a new variable, the temperature T, is introduced as

$$\epsilon +p=\sigma T+\sum _{i=1}^j{\mu }_i{n}_i.$$

(8)

In what follows, we do not consider the effects of conserved charges and we assume that the corresponding chemical potentials vanish, ${\mu }_i=0$.

Let us note here that we intended to search for self-similar solutions, and in the next section we actually describe a rich family of new and exact solutions that have a trivial proper-time dependence but also depend on the coordinates predominantly through a certain scaling variable s. At the end, we describe solutions that break self-similarity, nevertheless they do it with an explicit factor that vanishes around mid-rapidity, so the concept of self-similarity and its actual breaking will make the investigation of the question of self-similarity of the solution an interesting one.

The definition of the scaling variable s is that it has a vanishing co-moving derivative:

$$u_{\mu }{\partial }^{\mu }s=0.$$

(9)

If s is a scaling variable that satisfies the above equation, then obviously any $s^{\prime }=s^{\prime }\left(s\right)$ that is a function of s only may also be introduced as a new scaling variable, given that

$$u_{\mu }{\partial }^{\mu }{s}^{\prime }={\partial }_s{s}^{\prime }{u}_{\mu }{\partial }^{\mu }s=0.$$

(10)

In 1$+$1 dimension, let us rewrite the equations of relativistic hydrodynamics with the help of the Rindler coordinates $(\tau ,{\eta }_x)$ for the temperature T and for the rapidity of the fluid $\mathrm{\Omega }$. The longitudinal proper-time, referred to in what follows simply as the proper-time ($\tau$), the coordinate-space rapidity ${\eta }_x$ and the fluid rapidity $\mathrm{\Omega }$ are defined as

$$\begin{array}{cc}\hfill \tau & =\sqrt{t^2-r_z^2},\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill {\eta }_x& =\frac{1}{2}\ln \left(\frac{t+r_z}{t-r_z}\right),\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \mathrm{\Omega }& =\frac{1}{2}\ln \left(\frac{1+v_z}{1-v_z}\right).\hfill \end{array}$$

(13)

The fluid rapidity $\mathrm{\Omega }$ relates to the four-velocity and to the three-velocity as

$$\begin{array}{ccc}\hfill u^{\mu }& =& \left(\cosh \left(\mathrm{\Omega }\right),\sinh \left(\mathrm{\Omega }\right)\right),\hfill \end{array}$$

(14)

$$\begin{array}{ccc}\hfill v_z& =& \tanh \left(\mathrm{\Omega }\right).\hfill \end{array}$$

(15)

For notational clarity, let us also highlight here the definition of the observables, which are determined by measuring particle tracks in momentum-space. The pseudorapidity ${\eta }_p$ and the rapidity y of a final state particle with mass m and four-momentum $p^{\mu }$ are defined as

$$\begin{array}{cc}\hfill {\eta }_p& =\frac{1}{2}\ln \left(\frac{p+p_z}{p-p_z}\right),\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill y& =\frac{1}{2}\ln \left(\frac{E+p_z}{E-p_z}\right),\hfill \end{array}$$

(17)

where the modulus of the three-momentum is denoted by $p=\left|\mathbf{p}\right|=\sqrt{p_x^2+p_y^2+p_z^2}$.

In the temperature independent speed of sound or $\kappa \left(T\right)=const$ approximation, the relation between the pressure, temperature and entropy density simplifies as follows:

$$p=\frac{T\sigma }{1+\kappa }.$$

(18)

From now, let us assume that the fluid rapidity depends only on the coordinate space rapidity, but not on the proper time, so that

$$\begin{array}{ccc}\hfill \mathrm{\Omega }& =& \mathrm{\Omega }\left({\eta }_x\right).\hfill \end{array}$$

(19)

With the help of the above assumption, the relativistic energy and Euler equations can be rewritten for the temperature and for the fluid rapidity in terms of the Rindler coordinates, using the notation $T=T(\tau ,{\eta }_x)$ and $\mathrm{\Omega }=\mathrm{\Omega }\left({\eta }_x\right)$:

$$\begin{array}{ccc}\hfill {\partial }_{{\eta }_x}\mathrm{\Omega }+\kappa \left(\tau {\partial }_{\tau }+tanh\left(\mathrm{\Omega }-{\eta }_x\right){\partial }_{{\eta }_x}\right)ln\left(T\right)& =& 0,\hfill \end{array}$$

(20)

$$\begin{array}{ccc}\hfill {\partial }_{{\eta }_x}ln\left(T\right)+tanh\left(\mathrm{\Omega }-{\eta }_x\right)\left(\tau {\partial }_{\tau }ln\left(T\right)+{\partial }_{{\eta }_x}\mathrm{\Omega }\right)& =& 0.\hfill \end{array}$$

(21)

In this way, we obtained a set of partial differential equations for the temperature and for the fluid-rapidity using the variables $(\tau ,{\eta }_x)$. The solutions of these hydrodynamical equations are presented in the next section.

Let us clarify that we do not present the details of the derivation of these solutions in this manuscript: these complicated details go well beyond the scope of the presentation of the results. However, the validity of these solutions is straigthtforward to check by substituting them to the temperature and the Euler equations of Equations (20) and (21) and to the entropy equation of Equation (2). Note also that the above temperature and Euler equations, Equations (20) and (21) were derived in Rindler coordinates before, in arbitrary d number dimensions, as given in Appendix A of Ref. [60].

4. A New Family of Exact Solutions of Relativistic Hydrodynamics

We have discovered the following family of exact and explicit solutions of the equations of relativistic hydrodynamics:

$$\begin{array}{ccc}\hfill {\eta }_x\left(H\right)& =& \mathrm{\Omega }\left(H\right)-H,\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill \mathrm{\Omega }\left(H\right)& =& \frac{\lambda }{\sqrt{\lambda -1}\sqrt{\kappa -\lambda }}\arctan \left(\sqrt{\frac{\kappa -\lambda }{\lambda -1}}\tanh \left(H\right)\right),\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill \sigma (\tau ,H)& =& {\sigma }_0{\left(\frac{{\tau }_0}{\tau }\right)}^{\lambda }{\mathcal{V}}_{\sigma }\left(s\right){\left[1+\frac{\kappa -1}{\lambda -1}{\sinh }^2\left(H\right)\right]}^{-\frac{\lambda }{2}},\hfill \end{array}$$

(24)

$$\begin{array}{ccc}\hfill T(\tau ,H)& =& T_0{\left(\frac{{\tau }_0}{\tau }\right)}^{\frac{\lambda }{\kappa }}\mathcal{T}\left(s\right){\left[1+\frac{\kappa -1}{\lambda -1}{\sinh }^2\left(H\right)\right]}^{-\frac{\lambda }{2\kappa }},\hfill \end{array}$$

(25)

$$\begin{array}{ccc}\hfill \mathcal{T}\left(s\right)& =& \frac{1}{{\mathcal{V}}_{\sigma }\left(s\right)},\hfill \end{array}$$

(26)

$$\begin{array}{ccc}\hfill s(\tau ,H)& =& {\left(\frac{{\tau }_0}{\tau }\right)}^{\lambda -1}\sinh \left(H\right){\left[1+\frac{\kappa -1}{\lambda -1}{\sinh }^2\left(H\right)\right]}^{-\lambda /2}.\hfill \end{array}$$

(27)

A new property of the above equations is that the space-time rapidity ${\eta }_x$ dependence of the hydrodynamic fields are given as parametric curves with explicit proper-time dependence. If one intends to analyse these solutions at a given proper-time $\tau$, one may plot for example the $({\eta }_x\left(H\right),F(\tau ,H))$ parametric curve, where F stands for any of the hydrodynamic fields: $F=\left\{T,\sigma ,p,\epsilon ,\mathrm{\Omega },v_z,\dots \right\}$. The parameter of these curves, $H=\mathrm{\Omega }-{\eta }_x$, has a clear physical meaning; it is the difference of the fluid rapidity $\mathrm{\Omega }$ and the space-time rapidity ${\eta }_x$. Some special properties of these parametric curves is analyzed in the next subsection.

The quantity $\lambda$ is a constant of integration that is a measure of the relativistic acceleration, and as we show subsequently, this parameter also controls the width of the rapidity distribution and it can be determined from the fits to experimental data. It is worth noting here that these solutions are obtained in the physical region, where $1\le \lambda <\kappa$, given that experimental data indicate $\lambda \simeq$ 1.05–1.2 at RHIC and LHC energies, while in 200 GeV Au + Au collisions, the average value of the speed of sound was measured by the PHENIX collaboration, yielding $\kappa =1/c_s^2\approx 10$ [68]. The constants of normalization, ${\sigma }_0$ and $T_0$ are chosen to denote the initial conditions, the value of the entropy density and temperature at the initial proper-time $\tau ={\tau }_0$ and at mid-rapidity, where ${\eta }_x=\mathrm{\Omega }=H=0$. The scaling functions for the entropy and the temperature profile are denoted as ${\mathcal{V}}_{\sigma }\left(s\right)$ and $\mathcal{T}\left(s\right)$, but, as indicated in the above solution, only one of them (for example, the temperature profile function) can be chosen arbitrarily. The scaling variable s above is normalized so that its proper-time dependent prefactor is unity at the initial proper-time, $\tau ={\tau }_0$ and its value vanishes at midrapidity, $s(\eta =0)=0$ as usual.

4.1. Discussion and Limiting Cases

One of the most important properties of the above equations is that they describe a family of exact solutions: for every positive definite, univariate scaling function of the temperature $\mathcal{T}\left(s\right)$ one finds a corresponding solution, so that (apart from an overall normalization factor) the temperature and the entropy density depends predominantly on the space-time rapidity ${\eta }_x$ through this scaling function $\mathcal{T}\left(s\right)$. This property is inherited from the 1 + 1 dimensional family of solutions of relativistic hydrodynamics described by Csörgő, Grassi, Hama and Kodama in Refs. [36,38], that are recovered in the $H\ll 1$ and $\lambda \to 1$ limit of our new solutions. However, in our case, additional factors are also present in the solutions that break their self-similarity explicitly.

The space-time rapidity dependence of our new solutions is given by Equations (23)–(27). These equations describe parametric curves. This point requires a careful analysis as the physical solutions are to be given not as parametric curves, but as functions of the space-time rapidity. One way to handle this nature of the parametric curves is to limit their domain of applicability to those regions, where they correspond to functions. These regions are limited by the points, where the derivative of the space-time rapidity with respect to the parameter of the curve vanishes. For to this reason, we discuss the solutions only in a finite domain around mid-rapidity in this manuscript. From the criterion ${\partial }_H{\eta }_x\left(H\right)=0$, we get the following lower and upper limits for the applicability of our solution:

$$\begin{array}{cc}\hfill {\eta }_{min}& =-\frac{\lambda }{\sqrt{\left(\lambda -1\right)\left(\kappa -\lambda \right)}}\arctan \left(\sqrt{\frac{\kappa -\lambda }{\kappa \lambda -\kappa }}\right)+\mathrm{atanh}\left(\frac{1}{\sqrt{\kappa }}\right),\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill {\eta }_{max}& =\frac{\lambda }{\sqrt{\left(\lambda -1\right)\left(\kappa -\lambda \right)}}\arctan \left(\sqrt{\frac{\kappa -\lambda }{\kappa \lambda -\kappa }}\right)-\mathrm{atanh}\left(\frac{1}{\sqrt{\kappa }}\right),\hfill \end{array}$$

(29)

where ${\eta }_{max}=-{\eta }_{min}$. Thus, the domain of validity of our solutions is limited to the ${\eta }_{min}<{\eta }_x<{\eta }_{max}$ interval. However, the parametric curves are defined outside this interval, too. This suggests the possibility of the extension of our solutions, with the help of shock-wave equations, to the forward and backward rapidity regions. Their discussion is straightforward, but goes beyond the scope of the present manuscript.

The finiteness of the mid-rapidity region implies that there is an upper bound for the modulus of the possible values of the fluid rapidities:

$$\left|\mathrm{\Omega }\right|\le \frac{\lambda }{\sqrt{\lambda -1}\sqrt{\kappa -\lambda }}\arctan \left(\sqrt{\frac{\kappa -\lambda }{\kappa \left(\lambda -1\right)}}\right),$$

(30)

which implies that this solution is finite in terms of the fluid rapidity $\mathrm{\Omega }$. This also implies that the solution yields a finite observable rapidity and pseudorapidity distribution.

Let us also clarify that our solutions are in fact not self-similar: although they contain factors that depend on the space-time rapidity ${\eta }_x$ only through the scaling variable s, the solutions for the temperature and the entropy density contain additional factors that depend on ${\eta }_x$ through the parameter H as well. This implies that our solutions are approximately self-similar only in a small domain near mid-rapidity, but their self-similarity is explicitly violated by terms that become more and more important with increasing deviations from mid-rapidity.

Let us note that, in the $c_s^2=1/\kappa =1$ case, the solution described by Equations (23)–(27) reproduces the 1 + 1 dimensional CNC solution of Refs. [59,60], as expected. If $\kappa =1$ or if $H\ll 1$, the fluid rapidity is an explicit function of the coordinate rapidity, given as

$$\mathrm{\Omega }\simeq \lambda {\eta }_x,\mathrm{for}{\eta }_x\ll 1/(\lambda -1).$$

(31)

However, let us also emphasize at this point that the $\kappa \to 1$ and the $\lambda \to 1$ limits are not interchangeable: The CNC limit corresponds to the $\kappa \to 1$ limit of the extension of the above solutions to the $1\le \kappa <\lambda$ domain of the model parameters, however the detailed description of this class of solutions is beyond the scope of the present manuscript. In the $\kappa =1$ limit, we recover the CNC, while in the $\lambda =1$ limit, we recover the CGHK family of solutions.

Let us also clarify, that our new family of solutions includes the Hwa–Bjorken asymptotic solution as well. To see this, one has to proceed carefully, as the $H\to 0$ and the $\lambda \to 1$ limits of our new family of solutions are not interchangeable. The Hwa–Bjorken solution is $\mathrm{\Omega }={\eta }_x$, and it corresponds to the $H=0$ and $\lambda \to 1$ from above limiting case. As $H\to 0$, $\mathrm{\Omega }\to \lambda {\eta }_x$, and if we take the $\lambda \to 1$ limit after the $H\to 0$ limit, we recover the flow field of the Hwa–Bjorken solution. In this case, the domain of validity, the interval $({\eta }_{min}\le {\eta }_x\le {\eta }_{max})$ also approaches the whole space-time rapidity region, $(-\infty <{\eta }_x<\infty )$. After these steps, we can consider the $\mathcal{T}\left(s\right)=1$ special case, to demonstrate that our new family of solutions includes the Hwa–Bjorken solution as a special limiting case. However, if one takes the $\lambda \to 1$ limit too early, before the $H\to 0$ limit, one could find an unphysical, divergent hence uninteresting limiting case.

4.2. Graphical Illustrations

The temperature map of our new family of solutions is indicated in Figure 1. The finiteness of these solutions means that they are defined in a cone, $({\eta }_{min}\le {\eta }_x\le {\eta }_{max})$, that lies within the forward light cone. The pseudorapidity dependence of the temperature at various constant values of the proper-time is indicated on Figure 2. It is remarkable, how similar these exact solutions are to the figures of those Monte-Carlo simulations that mapped out the space-time picture of heavy ion collisions, see for one of the first examples the figures of Refs. [69,70].

The map of the fluid-rapidity distribution of our new family of solutions, $\mathrm{\Omega }(t,r_z)$, is illustrated in Figure 3. This plot shows not only the finiteness of these solutions, which corresponds to a cone within the forward light cone, but also that in our case the fluid rapidity is a function of the space-time rapidity only, so that $\mathrm{\Omega }=\mathrm{\Omega }\left({\eta }_x\right)$, independently of the proper-time $\tau$. The space-time rapidity dependence of the fluid-rapidity at various constant values of the proper-time is indicated on Figure 4. This figure also indicates that our new exact solutions, although formally different, numerically are very close to the $\mathrm{\Omega }=\lambda {\eta }_x$ flow field, that is to the CNC solutions of Refs. [59,60].

In these earlier CNC solutions, the $\mathrm{\Omega }=\lambda {\eta }_x$ flow rapidity distribution was an exact solution, but only for a super-hard, hence unrealistic equation of state, corresponding to the $\kappa =1$, or speed of sound $c_s=1$ equation of state. Although this equation of state has been changed drastically, the speed of sound was decreased from speed of light to any temperature independent value, the flow field changed surprisingly little: most of its change is only at forward and backward rapidities and it is only a few percent for the experimentally relevant values of the parameters $\lambda$ and $\kappa$.

As the acceleration parameter $\lambda \to 1$, the space-time rapidity region opens up, as it approaches the whole horizontal axes. At the same time, the temperature approaches a constant value that becomes independent of the value of the space-time rapidity coordinate ${\eta }_x$: we recover the boost-invariant Hwa–Bjorken solution. This is a common property of Figure 1, Figure 2, Figure 3 and Figure 4 .

5. Observables: Rapidity and Pseudorapidity Distributions

To be able to compare our solution to data, let us calculate analytically the observable quantities, that correspond to the rapidity and to the pseudorapidity distributions, both defined in momentum space. The detailed calculation is not included in this manuscript, but the main steps are summarized as follows. We assume the simplest possible form of the temperature scaling function, using $\mathcal{T}\left(s\right)=1$. We utilize a Boltzmann approximation to evaluate the phase-space densities on the freeze-out hypersurface, which is assumed to be pseudo-orthogonal to the four-velocity $u^{\mu }$. Consequently, the normal vector of the freeze-out hypersurface is parallel to $u^{\mu }$, i.e., $d{\sigma }^{\mu }=u^{\mu }d{\eta }_x$, where $d{\sigma }^{\mu }$ is the infinitesimal form of the normal vector of the freeze-out hypersurface:

$$d{\sigma }^{\mu }=\frac{1}{A\left({\eta }_x\right)}\left({\partial }_{{\eta }_x}{r}_z,{\partial }_{{\eta }_x}t\right)d{\eta }_x,$$

(32)

where $A\left({\eta }_x\right)$ is a normalization factor [71]. We find that the equation of the freeze-out hypersurface is the following:

$$\frac{\tau \left(H\right)}{{\tau }_f}={cosh}^{\frac{\kappa }{\lambda -\kappa }}\left(H\right){\left[1+\frac{\kappa -1}{\lambda -1}{sinh}^2\left(H\right)\right]}^{\frac{\lambda }{2(\kappa -\lambda )}}.$$

(33)

Here, ${\tau }_f$ stands for the proper-time of the kinetic freeze-out at midrapidity, corresponding to the $H={\eta }_x=\mathrm{\Omega }=0$ parameter values. During the calculation of the rapidity distribution, we assume that we are not too far away from the midrapidity region, and with the help of a saddle-point approximation we perform the integral over ${\eta }_x$ to obtain the following formula:

$$\frac{dn}{dy}\approx {\left.\frac{dn}{dy}\right|}_{y=0}{cosh}^{-\frac{1}{2}\alpha \left(\kappa \right)-1}\left(\frac{y}{\alpha \left(1\right)}\right)exp\left(-\frac{m}{T_f}\left[{cosh}^{\alpha \left(\kappa \right)}\left(\frac{y}{\alpha \left(1\right)}\right)-1\right]\right),$$

(34)

where y is the rapidity, and $\alpha \left(\kappa \right)$ stands for

$$\alpha \left(\kappa \right)=\frac{2\lambda -\kappa }{\lambda -\kappa }.$$

(35)

A similar, but $\kappa$ independent exponent, $\alpha (\kappa =1)\equiv \alpha \left(1\right)$ was denoted with $\alpha$ in the CNC solutions of Refs. [59,60].

The normalization factor of this distribution is

$${\left.\frac{dn}{dy}\right|}_{y=0}=\frac{R^2\pi {\tau }_f}{{\left(2\pi \hslash \right)}^3}\sqrt{\frac{{\left(2\pi T_{f}m\right)}^3}{\lambda (2\lambda -1)}}exp\left(-\frac{m}{T_f}\right),$$

(36)

where the $R^2\pi$ multiplication factor is the transverse area and m is the mass of the observed particle (predominantly pion).

We emphasize that the result of Equation (34) reproduces the flat Hwa–Bjorken rapidity plateaux [17,18] in the $\lambda \to 1$ from above limit. In the $\kappa \to 1$ limit, our calculations also reproduce the rapidity distribution of the CNC solution of Refs. [59,60]. Using a mean value theorem for definite integrals and a saddle point approximation, the pseudorapidity distribution was found to be proportional to Equation (34), and the constant of proportionality to be the average value of the $\mathcal{J}=\frac{dy}{d{\eta }_p}$ Jacobi determinant. In this way, the pseudorapidity distribution can be expressed as a parametric curve $({\eta }_p\left(y\right),\frac{dn}{d{\eta }_p}\left(y\right))$ where the parameter of the curve is the rapidity y, and the pseudorapidity distribution is given by the following form:

$$\frac{dn}{d{\eta }_p}\left(y\right)\approx {\left.\frac{dn}{dy}\right|}_{y=0}\frac{\langle p_T\left(y\right)\rangle cosh\left({\eta }_p\right)}{\sqrt{m^2+{\langle p_T\left(y\right)\rangle }^{2}cosh\left({\eta }_p\right)}}{cosh}^{-\frac{1}{2}\alpha \left(\kappa \right)-1}\left(\frac{y}{\alpha \left(1\right)}\right)exp\left(-\frac{m}{T_f}\left[{cosh}^{\alpha \left(\kappa \right)}\left(\frac{y}{\alpha \left(1\right)}\right)-1\right]\right),$$

(37)

Note that ${\eta }_p\left(y\right)$ is the pseudorapidity evaluated as a function of rapidity y with the help of the rapidity dependent average transverse momentum, $\langle p_T\left(y\right)\rangle$ defined as:

$$\langle p_T\left(y\right)\rangle \approx \frac{\sqrt{T_f^2+2m{T}_f}}{{\displaystyle \phantom{(}1+\frac{{\displaystyle \phantom{(}\alpha \left(\kappa \right)}}{{\displaystyle \phantom{(}2\alpha {\left(1\right)}^2}}\frac{{\displaystyle \phantom{(}{T}_f+m}}{{\displaystyle \phantom{(}{T}_f+2m}}{y}^2}}.$$

(38)

Note, that the same functional form, a Lorentzian shape was assumed for the rapidity dependence of the slope of the transverse momentum spectrum in the Buda–Lund hydro model of Ref. [41]. This form was obtained based on arguments of the analyticity of the inverse temperature profile in rapidity and the coefficient of the $y^2$ dependence was considered also very recently in Refs. [66,67] as a free fit parameter. Our exact solutions of fireball hydrodynamics express this coefficient with two observables, the freeze-out temperature $T_f$ and the mass, and the speed of sound parameter of the solution, $\kappa =1/c_s^2$. Thus, from the measurement of the rapidity dependence of the transverse momentum spectra, the equation of state parameter can in fact be fitted to data, at least in principle. Note that such an approximately Lorentzian rapidity dependence of the slope of the transverse momentum spectra has also been experimentally tested in h + p reactions at $\sqrt{s}=22$ GeV by the EHS/NA22 collaboration in Ref. [42], so apparently this feature of our exact solutions is realistic: it has experimental support.

The normalized rapidity and pseudorapidity distributions are illustrated in Figure 5. They look rather realistic and work has just been started to fit them to experimental data in proton–proton and heavy ion collisions at BNL RHIC and at CERN LHC energies. Some of the first results are shown in the next section.

6. First Comparisons to p + p Data at RHIC and LHC

The first comparisons of our new, exact solutions of relativistic hydrodynamics to experimental data are shown in Figure 6 and Figure 7.

Apparently, CMS data on $dn/d{\eta }_p$ pseudorapidity densities, as measured in $\sqrt{s}=$ 7 [72] and 8 TeV [73] p + p collisions at LHC, are described with Equations (37) and (38), as obtained from our exact family of solutions, corresponding to Equations (23) and (27). Let us also clarify that the two investigated CMS datasets were measured with somewhat different trigger conditions at $\sqrt{s}=7$ and 8 TeV, which explains the difference between their absolute normalizations.

The confidence level of these fits, CL > 0.1 % indicates that our new solution is not inconsistent with these CMS data, and the parameters can be interpreted as meaningful ones. However, given the large errors of the measurement, the parameters of the $\sqrt{s}=7$ and $\sqrt{s}=8$ TeV reactions cannot be distinguished based on data fitting, given that their values are, within three standard deviations and within the systematic error of this analysis, the same.

One should also be careful not to draw too strong conclusions here because our solutions are limited in rapidity range. The domain of the validity of our solution is given by Equations (28) and (29) and depending on the values of the fit parameters $\kappa$ and $\lambda$ these range from $|{\eta }_x|<1.03$ to $|{\eta }_x|<2.5$. Given that our solutions are finite, the propagation of shock-waves should be investigated in the target and in the projectile rapidity region, but instead we limit the range and plot the comparison to data only in the limited $|{\eta }_p|<2.5$ pseudorapidity interval.

Let us also note that we can describe the data with physically very different equations of state. For example, using the equation of state $\epsilon =\kappa p$ with a realistic, $\kappa =10$ value, in agreement with $c_s^2=1/\kappa \approx 0.1$ of Ref. [68], we get as good a description in Figure 6 as with the ideal, massless relativistic gas equation of state, corresponding to $\kappa =3$ in Figure 7. Similarly good quality fits were obtained by members of our group to describe the same pseudorapidity distribution, but measured in the complete pseudorapidity range, using the unrealistic, superhard equation of state $\kappa =1$ [66]. Thus, one of the apparent conclusions is that the pseudorapidity distributions are not very sensitive, at the present level of their precision, to the equation of state parameter $\kappa$ as long as the flow field remains approximately linear, $\mathrm{\Omega }\approx \lambda {\eta }_x$.

7. Summary

We have found a new family of analytic and accelerating, exact and finite solutions of perfect fluid hydrodynamics for 1$+$1 dimensionally expanding fireballs. These solutions are defined in a cone within the forward lightcone, and can be used to fit rapidity or pseudorapidity densities at RHIC and at LHC. Near mid-rapidity, the resulting rapidity and pseudorapidity distributions from the new solution approximate well those calculations obtained from the earlier CNC solution, and near mid-rapidity, the longitudinal flow approximates the $v=tanh\left(\lambda {\eta }_x\right)$ CNC flow profile. Due to small correction terms, our solutions in principle allow determining both the acceleration parameter $\lambda$ and the equation of state parameter $\kappa$ directly from fitting the measured (pseudo)rapidity densities as well as the rapidity dependence of the mean transverse momentum high energy proton–proton and heavy ion collisions at RHIC and at LHC. The Hwa–Bjorken [17,18], CGHK [36,37] and CNC solutions [59,60] are recovered in appropriate limits of the model parameters.

Further generalizations of these solutions to 1$+$3 dimensional, transversally expanding and/or rotating exact solutions, as well as to solutions with a temperature dependent speed of sound, are desirable and are being explored at the time of closing this manuscript.