Family of New Exact Solutions for Longitudinally Expanding Ideal Fluids

We report on the discovery of new analytical solutions of the equations of relativistic ideal hydrodynamics. In this solution, the fluid expands in the longitudinal direction and contains a plateau structure that extends over a finite range in rapidity and can be either symmetric or asymmetric in that variable. We further calculate the corresponding pseudo-rapidity distribution of hadron yields, and find decent agreement with experimental measurements in high-energy Pb+Pb, Au+Au, p+Pb, and d+Au collisions.

Introduction -Relativistic heavy-ion collisions allow systematic laboratory-based studies of a color-deconfined phase of matter -the Quark-Gluon Plasma(QGP). Owing much to the vigorous program pursued at the Relativistic Heavy-Ion Collider (RHIC) and at the Large Hadron Collider (LHC), one of the breakthroughs in theoretical relativistic heavy-ion physics has been the realization of the great success of numerical hydrodynamic simulations in describing the evolution of QGP, as well as understanding and predicting experimental measurements highlighting the collective behavior of the observed hadrons (see e.g. Refs. [1][2][3]). Experimentally, this collective behavior is observed through measurements and analyses of multi-particle correlation. On the theory side, fluid dynamics governs the time evolution of T µν , the energy-momentum tensor. More specifically, -and anticipating the use of curvilinear coordinates -T µν evolves following the conservation laws where D µ is a covariant derivative and Γ µ ρµ the Christoffel symbol.
In parallel with the remarkable progress made in numerical fluid dynamics, the study of analytical solutions remains useful in capturing intuitive pictures and important features. In that context, Landau, Khalatnikov, and Belenkij gave the first implicit solution formulated for these equations in [4][5][6]. Later on, a simple solution was found independently by Hwa [7] and Bjorken [8], the latter formulation is now known as the Bjorken flow. The Bjorken flow depends only on the proper time τ , and it is invariant under a Lorentz boost along the expansion (longitudinal) direction. In recent decades, a new family of solutions for a longitudinal expanding fluid was found by Csörgő, Nagy, and Csanád [9], where the rapidity profile is symmetric. Other analytical formulation for the fluid dynamics of longitudinally expanding fluid now also exist [10][11][12][13][14][15].
There also have been some developments in finding analytical solutions with non-trivial transverse structure. * shuzhe.shi@stonybrook.edu Notably, taking the conformal Equation of State (EoS), a solution was found by Gubser [16] which is boostinvariant in the longitudinal direction and expands in the transverse plane. Another solution based on spherical expansion which allows non-trivial acceleration and rotation was found by Nagy [17], and more solutions with viscous effect were highlighted by Hatta, Noronha, and Xiao [18].
The QGP system evolving in relativistic heavy ion collisions is of course not boost-invariant. There is a finite range in rapidity that contains the hot medium, while the system is dilute outside this rapidity window. In this work, we introduce a family of new solutions to the 1+1D hydrodynamic equations, which is not boostinvariant and also can be either symmetric or asymmetric in rapidity. Starting from such a solution, we further compute the corresponding pseudo-rapidity distribution of hadron multiplicity frozen-out from the isothermal hypersurface. By choosing appropriate parameters, we find the pseudo-rapidity distribution computed from the analytic solution agrees reasonably well with experimental measurements [19][20][21][22][23].
Hydrodynamics in 1+1 D -We adopt the Milne coordinate system which combines the time and longitudinal coordinates into proper time, τ ≡ √ t 2 − z 2 , and spatial rapidity, η ≡ 1 2 ln t+z t−z . We focus on systems that are homogeneous in the transverse plane but contain nontrivial rapidity structure, which yields u x = u y = 0, and u η = 0. The relevant hydrodynamic equations are then: Furthermore, we neglect viscous corrections in the stress tensor so that T µν = (ε + p)u µ u ν − p g µν is the ideal fluid stress-energy tensor, and employ a simple equation of state p = c 2 s ε. Combining Eqs. (2) and (3) in two independent ways, one obtains the following two equations where ε 0 is a constant parameter with units , we introduce the fluid-rapidity ξ to express the velocity vector u τ = cosh ξ 2 and τ u η = sinh ξ 2 , so that the normalization condition, u µ u µ = 1, is automatically satisfied. With those new variables, Eqs. (4-5) may be expressed as Applying x − ∂ − to (6), and similarly x + ∂ + to (7) and subtracting the results, one can cancel out the ε-dependent terms and obtain the evolution equation only for ξ So far, no assumptions have been made. Eq. (8) and one of Eqs. (6) or (7) form a complete set of hydro equations that is equivalent to that of (2-3). From now on, we focus on the special case where the fluid-rapidity (ξ) can be separated as the superposition of an x + -dependent part and an x − -dependent part, i.e., With this ansatz, Eq. (8) can be simplified to Note that the left-hand-side of (10) is independent of x + , whereas the right-hand-side is independent of x − . The equality (10) can be fulfilled if and only if both sides equal to a constant, and we denote such a constant as β.
Combining the Eqs. (9-10) with the Eqs. (6-7), one finds For a general β, there is an analytic solution for ξ ± but not for ε. However, there exists a simple analytic solution in the case where β = 0, which is equivalent to the condition ∂ + ∂ − ln ε ε0 = 0, namely the energy density can be separated as the production of x + -and x − -dependent parts, So far, we have simplified the equations to be solved by focusing on the systems where fluid-rapidity and energy density can be separated into x + -and x − -dependent parts (9,12), -the applicability to heavy-ion collisions will be justified later on -and we found the following family of solutions for the flow rapidity where η 0 and a are dimensionless constants, t 0 has the unit of time. Substituting Eq. (13) into Eqs. (6-7) and solving the resulting equations, we obtain where ε 0 has been re-defined to absorb extra constants. Finally, we express the above solutions, i.e. the energy density and velocity field, in Milne coordinates: Here, we discuss the meaning of the parameters appearing in the solution: • The parameter η 0 simply shifts the spatial rapidity, and can always be absorbed by applying a Lorentz transformation τ → τ , η → η + η 0 . Hence, one can set η 0 = 0 with no loss of generality.
• τ 0 is a positive constant that scales the proper time τ , but is not necessarily its initial value. In other words, the hydro evolution can start from τ /τ 0 < 1.
• t 0 is a non-negative constant with units of time. It controls the width of the rapidity structure and can be treated as the starting time of the Bjorken-like expansion. See below for explanation.
• a quantifies the asymmetry in rapidity of the solution. It covers the range In particular, the solutions corresponding to a = A (where A is some arbitrary value) and a = 1/A are parity reflections (η ↔ −η) of each other (see Fig. 1 a,b). When a = 1, the solution is symmetric (see Fig. 1 c), and the energy density takes a simple form This solution is a generalization of the Bjorken flow which re-defines the time t → t + t 0 . In doing so, we obtain a rapidity-dependent solution. In some sense, Eq. (18) can be regarded a Bjorken-like solution when the beam energy is finite and hence, the overlap time (t 0 ) is finite. The solution (15)(16)(17) returns to the Bjorken solution when a = 1 and t 0 = 0. In the other extreme limit when a = has the appearance of a sigmoid, or a smooth step function (Fig. 1 (d)). To obtain the solution Eqs. (15)(16)(17), the only assumption that has been made was that the fluid-rapidity and the energy density can be separated into x + -and x −dependent parts, i.e. Eqs. (9,12). Noting that x + = τ e η / √ 2 corresponds to the target(backward) contribution, whereas x − = τ e −η / √ 2 corresponds to the projectile(forward) contribution, one can interpret Eqs. (9,12) as the separation of target and projectile contributions: Here the subscript P denotes the projectile and the sub-script T denotes the target. We point out that Eq. (20) can be realized in Glasma-based models: in Refs. [24,25], it was argued that the energy deposition in the transverse plane, averaged over color-charge fluctuations, can be separated as the product of saturation scales of the target and the projectile: Hence, one may interpret the x ± contributions to the energy density in Eq. (20) as the target and projectiles saturation scales: For values of rapidity such that e ∓η t 0 /τ , the above expressions approaches the solution of the JIMWLK evolution [26] with a constant speed (λ ∓ ≡ A similar rapidity dependence is also obtained in Ref. [27], which is based on parton saturation and classical Chromo-Dynamics. The asymmetry parameter a = 1 implies different evolution speed for projectile and target in asymmetric collisions. Also note that a positive t 0 prevents the appearance of a divergence at large η, i.e. near the source nuclei. In Fig. 2, we separately plot the target and projectile contributions to both the fluid rapidity and the energy density. The solution, Eqs. (22)(23), exhibits a plateau in the energy density near the source nuclei, followed by an exponential tail at large distance in spatial rapidity. In an asymmetric system (a = 1.02), both the heights of the energy-plateau and the slopes of the exponential tail are different for the projectile and target.
Hadron Distribution -It would be interesting to examine the applicability of the solution (15)(16)(17) to the medium created in relativistic heavy-ion collisions. A natural criteria is the rapidity dependence of particle yield, which can be computed in the theory and also measured in experiments.
One can employ the Cooper-Frye formula [28] to calculate the momentum distribution of observed hadrons from a given hydrodynamic profile: with Σ f being the freeze-out hypersurface, d 3 σ µ the surface element, the +(−) sign is taken for baryons(mesons), and Θ(u·p) is a step function to ensure that particles always move out of the medium. Also, u · p = u µ p µ is the energy of the particle in the fluid cell rest frame. The hypersurface volume is given by (see Ref. [29] and references therein) where ζ, ζ , and ζ are the coordinates for the hypersurface.
After a tedious but straightforward calculation (see App. A), we obtain the pseudo-rapidity distribution of where S is a scaling factor proportional to the transverse area and can be adjusted according to the overall particle production rate, and T ini ≡ (ε 0 /ε f ) 1/4 T f indicates the initial temperature. The width (w) of the rapidity distribution and the slope at y p = 0 are found to be For details, see App. A. Taking the parameters τ ini = 0.4 fm/c, T f = 145 MeV, in line with phenomenological analyses [30,31] and setting t 0 = 0.1 fm/c in order to match the plateau width, we plot the resulting pseudorapidity distribution of charged multiplicity, dN ch /dη p , as the sum of π ± , K ± , and p(p) contributions as the solid curves in Fig. 3. For a satisfactory comparison with the experimental data [19][20][21][22][23], we set T ini /T f = 2.0, 2.0, 2.2, and 1.9 for p+Pb, d+Au, Pb+Pb, and Au+Au, respectively. Also, a = 1.05(1.00) for asymmetric(symmetric) systems. We have adjusted the overall scaling factor, S, and shifted the pseudorapidity -through the redefinition of reference frame -by −0.3 and −0.4 unit in p+Pb and d+Au comparison, respectively. We observe in Fig. 3 that the solid (theory) curves share compelling qualitative features with experimental results, although quantitative differences do remain. We note that Eq. (28) is the Cooper-Frye distribution for particles created at the freeze-out hypersurface, and remark that the discrepancy with data might be due to the absence of resonance decay and hadron scattering effects. Noting that both of these effects would smear out the rapidity distribution, we perform a rough estimation for the post-hadron-cascade distribution by convoluting the Cooper-Frye distribution (28) with a Gaussian smearing kernel, dN ] dN dηp , with a width of σ = 1. The smeared distributions are shown as dashed curves in Fig. 3, which now exhibit reasonable agreement with the experimental data. While it is understood that some important features of the collision dynamics are not treated here -like the pre-equilibrium phase and viscous behavior -it is revealing that overall features of the medium created in symmetric and asymmetric heavy-ion collisions can be reproduced by our exact hydrodynamic solution. Interestingly, it may well turn out that disagreement with data is more interesting than agreement in this case, and may be used to signal departure from ideal fluid-dynamical behavior.
Summary and Discussion -To summarize, we derived a new family of exact solutions to the 1+1D ideal hydrodynamics that can describe heavy ion collisions at finite collision energies. A solution can be either symmetric or asymmetric, and it is contained within a finite rapidity range. Based on such a solution, we further computed the distribution of final state hadrons. By taking appropriate value for the parameters in the solution, we found reasonably good agreement with the experimental measurements in relativistic d+Au, p+Pb, Au+Au, and Pb+Pb collisions. With its key property summarized in Eqs. (29)(30), this flexible solution should be useful in providing guidance for phenomenological modeling of the longitudinal initial condition of heavy-ion collisions. In addition, exact solutions are valuable in the calibration of numerical integrations of hydrodynamical equations [32].
Moreover, the "discovery" of generalized Bjorken flow, Eq. (18), inspires us to point out a way to generalized any given solution of hydrodynamic equations. We note that hydrodynamic equations are covariant under translation in Minkowski spacetime coordinates -suppose T µν (x) satisfies the conservation equation ∂ µ T µν (x) = 0, then ∂ µ T µν (x + x 0 ) = 0 is valid for any constant x 0 . Therefore, given any boost-invariant solution, one can always perform a time translation by t → t = t + t 0 in the Minkowski coordinate, and the translated profile has non-trivial rapidity dependence and is also a solution to the hydrodynamic equations. While this is straightforward from a mathematical point-of-view, we emphasize that it leads to non-trivial physical consequence, particularly for the application in heavy-ion collisions. In the hydrodynamic simulation of heavy-ion collisions, one typically assumes the whole system is initialized at given constant proper time, and the translation in Minkowski time leads to the following transformation in Milne coordinates, τ = (τ 2 +2t 0 τ cosh η+t 2 0 ) 1 2 . Hence, initialization at constant τ is equivalent to a different "initialization" scheme in τ and η, and leads to different hadron distributions.