Analytical insights into the interplay of momentum, multiplicity and the speed of sound in heavy-ion collisions

We introduce a minimal model of ultracentral heavy-ion collisions to study the relation between the speed of sound of the produced plasma and the final particles' energy and multiplicity. We discuss how the particles' multiplicity $N_{\textrm{tot}}$ and average energy $E_{\textrm{tot}}/N_{\textrm{tot}}$ is related to the speed of sound $c_s$ by $c_s^2=d \ln (E_{\textrm{tot}}/N_{\textrm{tot}})/d\ln N_{\textrm{tot}}$ if the fluid is inviscid, its speed of sound is constant and all final particles can be measured. We show that finite rapidity cuts on the particles' multiplicity $N$ and energy $E$ introduce corrections between $c_s^2$ and $d \ln (E/N)/d\ln N$ that depend on the system's lifetime. We study analytically these deviations with the Gubser hydrodynamic solution, finding that, for ultrarelativistic bosons, they scale as the ratio of the freezeout temperature $T_{\mathrm{FO}}$ over the maximum initial temperature of the fluid $T_{0}$; the non-thermodynamic aspect of these corrections is highlighted through their dependence on the system's initial conditions.


I. INTRODUCTION
Ultrarelativistic heavy-ion collisions can produce strongly-coupled quark-gluon plasma, a state of matter where the fundamental degrees of freedom of Quantum Chromodynamics (QCD) are deconfined [1][2][3][4].The plasma expands at velocities of the order of the speed of light, and its evolution can be described with relativistic hydrodynamic models [2,[5][6][7][8].Recombining into hadrons as it cools down, the plasma's properties can be studied from the resulting final particles that reach the detectors.
The primary connection point between the properties of the plasma and the equations of fluid dynamics is the equation of state.At zero or small values of baryon chemical potential µ B , the equation of state of nuclear matter has been computed ab initio using lattice QCD simulations [9][10][11][12][13].The zero-µ B lattice equation of state, matched at low temperature to a hadron resonance gas [3,9,10,[14][15][16], is almost universally used as a pre-determined input for hydrodynamic simulations of high-energy heavy-ion collisions. 1 While there are known subtleties in matching the lattice equation of state to a hadron resonance gas [20,21], these issues are likely much smaller than the broader uncertainties present in heavy-ion collision modeling.Consequently, there have been few recent attempts to extract the equation of state from heavy-ion measurements, with most efforts focusing on other properties of QGP like its viscosity. 2One notable extraction of the equation of state from heavy-ion measurements is ref.[22]; it compared heavy-ion simulations to measurements from the Relativistic Heavy-Ion Collider (RHIC) and the Large Hadron Collider (LHC) using Bayesian parameter inference, constraining simultaneously the equation of state and other model parameters, and finding results consistent with lattice calculations for the equation of state.
Building on earlier works [23][24][25], it was recently proposed in refs.[26,27] that the QGP equation of state can be measured more directly, by using the relation between the mean transverse momentum ⟨p T ⟩ and particle multiplicity dN ch /dη from ultra-central heavy-ion collisions to extract the speed of sound c s : where T eff is an effective temperature of an associated uniform fluid at rest possessing the same total energy and entropy as the QGP at freezeout.The argument behind eq.(1) in refs.[26,27] is that the speed of sound is defined as c 2 s ≡ dP/dε = d ln T /d ln s, the effective temperature is related to the mean transverse momentum of particles, T eff ≈ (1/3)⟨p T ⟩ [28], and that the entropy density is related to the particle multiplicity, s ∝ N ch .
Equation (1) has been employed by CMS collaboration [29] to extract results for the speed of sound that were found to be in very good agreement with lattice QCD calculations [30].Subsequently, ref. [31] employed numerical simulations of heavy-ion collisions to study mechanisms other than the equation of state that can influence the ratio d ln⟨p T ⟩/d ln [dN ch /dη], identifying initial state fluctuations, energy deposition ansätze and centrality selection as factors.In ref. [32], the authors of the original proposal assessed the robustness of eq. ( 1) with respect to changes in the equation of state, finding support in numerical simulations for eq.(1).
In this work, we study mathematically and numerically the validity of eq.(1) using inviscid relativistic fluid dynamics with a constant speed of sound.Because the speed of sound is assumed not to vary with temperature, there is no ambiguity on which value of the speed of sound should be recovered precisely, allowing for precision studies; in particular, our results are independent of the definition of the effective temperature T eff .Given that the arguments presented in refs.[26,27] do not involve any details of the system of interest, they should also be applicable in our simplified setting.Thus, we should be able to use our simplified model to test the essential conceptual content of the prescription.
We focus on initial temperature profiles that are fixed in shape, but whose overall normalization varies, akin to ultracentral heavy-ion collisions where the size of the energy deposition is fixed but the amount of deposited energy changes.Variations of the normalization of the initial temperature profile lead to changes in the final particles' multiplicity and average energy, which we show in Section II to be related to the speed of sound by where E tot is the total energy contained in the fluid and N tot is the total number of particles, while N mr = dN/dy| y=0 is the particle yield at midrapidity and ⟨p T ⟩ their midrapidity average transverse momentum.The first equality is demonstrated to be almost exact for inviscid fluids with a constant speed of sound.The second equality corresponds to the case where only a subset of final particles are measured -at midrapidity for example; in that case, the speed of sound is only recovered by taking the limit of asymptotically low freezeout temperature T FO .We corroborate this result in Section III with 3+1D relativistic hydrodynamic simulations, by scanning a range of rapidity cuts for the observables.
In Section IV, we use a simplified 1+1D relativistic hydrodynamic model (boost invariant and cylindrically symmetric [33,34]) to better understand the relation between the midrapidity multiplicity and average transverse momentum and the plasma's speed of sound.Consistent with eq. ( 2), we show numerically that this simple model can recover the speed of sound for sufficiently small  freezeout temperature T FO -fig.1a -but can also exhibit a strong dependence on the plasma's lifetime, as quantified by T FO , as shown in fig.1b. 3 Using the Gubser solution to inviscid relativistic hydrodynamics [35,36], we derive analytical expressions for the multiplicity and average transverse momentum of massless (or ultrarelativistic) bosons.From these results, we quantify explicitly the dependence on T FO of eq. ( 2) in the Gubser case, finding with T 0 the maximum initial state temperature of the plasma and q and τ 0 Gubser initial condition parameters.We discuss implications for heavy-ion collisions in the summary.Notation: We use units such that ℏ = c = k B = 1 and the mostly minus (+, −, −, −) metric signature.For a given arbitrary four-vector A µ , its cylindrical coordinate components are (A µ ) τ,r,ϕ,ηs = (A τ , A r , A ϕ , A ηs ).

II. TOTAL MULTIPLICITY, TOTAL ENERGY AND SPEED OF SOUND IN INVISCID HYDRODYNAMICS
The derivation of eq.(1) in refs.[26,27] is based on general thermodynamic arguments.In this section, we show that, for an inviscid fluid with constant speed of sound and specific variations of the initial conditions, it is possible to derive a nearly exact mathematical identity, whose form is close to eq. (1).We then discuss the mathematical conditions necessary for eq.(1) to hold when it is not possible to measure all final state particles.
We consider a relativistic inviscid fluid without conserved charges.The fundamental quantity describing the system is the energy-momentum tensor [37] T µν = εu µ u ν − P (g µν − u µ u ν ) , where ε is the energy density, u µ the flow velocity and P the pressure, which is related to the energy density ε by the equation of state P = P (ε).The dynamics of the system is given by the local conservation law which forms a closed system of partial differential equations for inviscid fluids.From eqs. ( 4) and ( 5) and the first law of thermodynamics, we also have the local conservation of the entropy current ∂ µ s µ = 0, where s µ = su µ , where T s = ε + P .

A. A one-parameter family of solutions
Let Ψ = {T (x ν ), u µ (x ν )} = {"temperature", "velocity"} be an arbitrary solution of the ideal fluid equations in 3+1 dimensions, with equation of state P = c 2 s ε, where c s is the (constant) speed of sound.Let us show that the one-parameter family of spacetime profiles Ψ(α) = {αT (x ν ), u µ (x ν )}, with α ∈ (0, +∞), are also solutions to the same ideal fluid equations.
First, let us express P and ε as functions of T .Differentiating P = c 2 s ε, we obtain sdT = c 2 s T ds, where s is the entropy density.This is a separable ordinary differential equation, which gives for some integration constant A. It follows that the entropy current and the stress-energy tensor are given by the following constitutive relations: Then, along the family Ψ(α) = {αT, u µ }, the entropy current and stress-energy tensor obey the following scaling law: Since α is a constant (in spacetime), we find that Considering that the hydrodynamic equations are just the conservation laws ∂ µ T µν = 0, we conclude that, if Ψ = {T, u ν } is a solution of ideal fluid dynamics, then the whole family Ψ(α) = {αT, u ν } is a set of equally acceptable solutions of the same equations.
B. An identity with the total energy and multiplicity Consider again the scaling law (8).If we integrate it along an arbitrary Cauchy surface Σ, we obtain where S tot and E tot are respectively the total entropy and the total energy (the latter evaluated in the laboratory frame).Thus, we immediately find that Note that, by Gauss' theorem, S tot is independent of the Cauchy surface Σ we chose (recall that the fluid is ideal, and thus ∂ µ s µ = 0).In particular, we can take Σ to be the freezeout surface Σ FO , upon which T = T FO .Then, considering that both the entropy density, s, and the particle number density, n, are functions only of the temperature, and thus are constant across Σ FO , we can take them out from the hypersurface integral (independently from the choice of equation of state).Hence, we have the following sequence of identities4 : where N tot is the total particle number at freezeout, which (according to Cooper-Frye [38]) coincides with the total multiplicity that is later detected. 5Note that the value of the ratio s(T FO )/n(T FO ) does not depend on α, because it is a pure function of T FO (which is a fixed parameter of the system, independent from α).Therefore, taking a derivative of the logarithm of the entropy eliminates the s(T FO )/n(T FO ) factor between S tot and N tot , yielding Thus, combining (10) with (12), we obtain the following identity: Note that this relation holds for arbitrary flow geometry in 3 + 1 dimensions: no symmetry assumptions are required.However, it must be kept in mind that the differential "d" is evaluated along the one-parameter family of states introduced in Section II A. None of the above relations hold for generic displacements.

C. From total observables to midrapidity observables
The identity (13) holds when the total energy and multiplicity can be measured.For the same result to hold with midrapidity measurements, it requires the midrapidity multiplicity and the transverse energy to share the same dependence on α as the corresponding total quantities, i.e.
where we remark that, for massless particles, the transverse energy can be expressed as dE/dy| y=0 = N mr ⟨p T ⟩ (see also Eq. ( 26)).Conditions ( 14) hold for a spherically symmetric system, for example, since particles and energy are emitted equally in all directions.Hence, the ratio between the number of transverse particles and the total number of particles is dictated by symmetry, and it does not depend on α.The same reasoning holds for the energy.Unfortunately, the general case (in the absence of symmetries) is not so straightforward.In the following, we will focus for clarity on the scaling behavior of the multiplicity, but a similar argument applies also to the energy.The quantity N mr , which is the number of particles emitted in the transversal direction, depends on α in two ways.First, when we rescale the solution by α, the particle current J µ is rescaled by a factor α c −2 s .This same effect enters also N tot , and these α-related factors cancel out.The second effect is that, along the solution Ψ(α), the freezeout surface is defined by the condition T (α) = αT (1) = T FO , which depends on α.This means, when we change α, we are sampling J µ on different surfaces.It follows that, instead of varying α, we may just set α = 1, and replace the derivative in α with a derivative in the freezeout temperature: Now, we note that the ratio N mr /N tot is only sensitive to the angular dependence of the flow, since it quantifies how many particles (with respect to the total) are traveling in the transverse direction.
It is reasonable to assume that, at very late times, when the shape of the flow has relaxed to its late-time profile, N mr (T FO )/N tot (T FO ) will approach a constant 6 .In other words, lim Applying the same reasoning to the transverse energy, we finally obtain the following claim: In the next sections, we shall confirm the ideas outlined in the present section employing a threedimensional hydrodynamic model of heavy-ion collisions and further explore the nonzero T FO effects with 1 + 1D models with cylindrical symmetry. 6There is a solid mathematical argument in support of this claim.Define φ = Nmr/Ntot, and view it as a function of the freezeout hypersurface.At sufficiently late times, when the transversal expansion takes over the longitudinal expansion, the function φ is monotonically non-decreasing.On the other hand, φ is bounded above by 1.A well-known theorem of analysis tells us that a bounded monotonic sequence must converge to a stationary value.Hence, φ approaches a positive constant as T FO approaches zero.

III. NUMERICAL RESULTS IN A SIMPLIFIED 3+1D MODEL OF ULTRARELATIVISTIC HEAVY-ION COLLISIONS
The equations of inviscid relativistic hydrodynamics can be written as by projecting eq. ( 5) along the flow velocity u µ (first equation) and orthogonal to it (second equation).We solve these equations numerically with MUSIC [39][40][41][42].The equation of state is again a constant speed of sound one: P = c 2 s ε.We use eq.( 18) in Milne coordinates τ = √ t 2 − z 2 , the hyperbolic time, and η s = tanh −1 (z/t) the spacetime rapidity.
We initialize hydrodynamics at τ = 0.4 fm.We use the optical Glauber model to collide two lead nuclei with zero impact parameter to generate the initial energy density profile on the transverse plane.The energy densities are then multiplied with a longitudinal profile function to obtain 3D initial conditions.We use the longitudinal profile from ref. [43], which has a flat region near midrapidity that tapers off following Gaussian tails.In this work, the flat region extends between ± 4 units in spatial rapidity η s .The width of the Gaussian tail is 0.6.The initial flow velocities in the transverse plane and in the η s direction are set to zero.
We convert the fluid's energy-momentum tensor T µν to particle degrees of freedom on a constant temperature hypersurface Σ FO at temperature T = T FO using with the Cooper-Frye formula [38,44] where g is the degeneracy factor, dΣ µ is the surface element of the constant temperature hypersurface and f (x, p) is the phase-space distribution of the particles.For simplicity, we consider that the particles generated are of a single species of massless bosons: We set the degeneracy factor to g = 1.
We calculate the slope observable d ln (E tot /N tot )/d ln N tot from the particle distribution obtained after the Cooper-Frye procedure.This is done for two different equations of state, both with constant speeds of sound: c 2 s = 1/3 and 1/4.We evaluate the slope observable at rapidity windows of size ∆y around midrapidity.Figures 2a and 2c show this ratio for the two different c 2 s .The calculation is done for two different p T cuts: we include all particles with p T < 2 GeV or with p T < 5 GeV.We see that for the largest rapidity window (∆y/2 = 10), when effectively all the particles are included, the true speed of sound is recovered.This statement is true regardless of the specific freezeout temperature chosen.This is the manifestation of the (almost) exact identity in eq. ( 13).The same identity can be verified in figs.2b and 2d, where the ratio of total energy per unit entropy escaping the freezeout hypersurface to the total entropy is given as a function of the size of spatial rapidity window ∆η s /2 around midrapidity.As we go to large enough spatial rapidity windows, all the energy and entropy is included and the true speed of sound is recovered.
The challenging issue is to recover the true speed of sound from the observables at midrapidity.From fig. 2, we see that we get closer to the true c 2 s as we go to lower and lower freezeout temperatures, which is consistent with eq. ( 17).At higher freezeout temperatures, the flow profile has not yet reached its late time limit, and the slope observable does not recover the correct c situation becomes worse when we open the rapidity window and take in more but not all particles from forward and backward rapidities.The initial temperatures are lower near the tail of the 3D distribution, causing them to freezeout even earlier, when they are farther from the late time flow profile.As the energy term is dominated by particles at large rapidities, the c 2 s estimation is further away from its true value.Only when all the particles are included is the extraction of the speed of sound correct, irrespective of the specific freezeout temperature.
Similarly, p T cuts affect the quality of c 2 s extraction at midrapidity, shown in figs.2a and 2c, as also pointed out in Ref. [32].Neglecting high-p T particles leads to an underestimate of the mean p T and hence underpredicts the speed of sound.In a typical ultracentral collision at RHIC or the LHC, T FO ≈ 100-150 MeV, and T 0 ≈ 400-600 MeV, implying T FO /T 0 ≈ 0.15-0.4; the results with T FO /T 0 ≥ 0.3 depicted in Fig. 2 show appreciable deviations between c 2 s and its estimator.

IV. NUMERICAL AND ANALYTICAL RESULTS IN A SIMPLIFIED 1+1D MODEL OF ULTRARELATIVISTIC HEAVY-ION COLLISIONS
If we focus on midrapidity measurements, we can approximate the fluid as invariant under boost along the collision axis of the nuclei, simplifying the equations of inviscid hydrodynamics (eq.( 18)).Since we are interested in ultracentral collisions, we make the further approximation that the fluid is radially symmetric in the plane transverse to the collision axis.With these two symmetries, the fluid four velocity can be parametrized as (u µ ) τ,r,ϕ,ηs = (cosh κ, sinh κ, 0, 0).In this setting, the two remaining independent hydrodynamic equations of motion can be cast as which can also be expressed in the simpler form [33,45] ∂ τ (T sinh κ) where we have employed the fact that for fluids with constant speed of sound, s = A(1+c 2 s ) T 1/c 2 s (see eq. ( 6)).Then, the numerical solutions of the above equations for the temperature and flow profiles, T (τ, r) and κ(τ, r), respectively, can be found.The Cooper-Frye particlization formula (19) can also be simplified in this symmetric setting [34]: where ϕ is the angle in the transverse plane, y is the particle's momentum rapidity, and m T = p 2 T + m 2 is the transverse mass for a particle of mass m, which reduces to m T = p T for the massless particles we study in this work.The integral Γ is a line integral on a curve Γ in the (τ, r)plane defined by T (τ, r) = T FO ; the curve Γ is the non-trivial part of the freezeout hypersurface, which must be determined from the solution of the hydrodynamic equations of motion.From the particle momentum distribution, the mean number of particles dN/dy and the mean particle momentum ⟨p T ⟩ at a given particle momentum rapidity y, are As shown in Appendix A, for massless particles produced at midrapidity, the formula for the multiplicity simplifies to while the transverse energy is given by where is the Euler gamma function [46].These formulas are used in the following sections to compute the observable (1), given solutions to inviscid hydrodynamic equations of motion, i.e. solutions for κ(τ, r) and T (τ, r).

A. Gaussian initial conditions
In order to compute the midrapidity particle yield N mr and the mean transverse momentum ⟨p T ⟩, we employ the following initial profiles for the temperature, We use σ = 5 fm as reference, based on the approximate radius of large atoms like lead and gold.We set the initial transverse velocity u r = sinh(κ) to a small value (10 −5 ) for numerical stability, effectively equivalent to using no initial transverse flow velocity.We solve the equations in Mathematica.In Fig. 3a, we first show the hypersurface curve Γ from eqs. ( 23), ( 25) and (26a).We see that, the lower the temperature, the more sharply peaked the curve is around r = τ when c 2 s = 1/3.It is from this region that most of the contribution to the integrals N mr and ⟨p T ⟩ originate.This will become clearer when we discuss Gubser flow in the next section.Evidently, for a fixed T 0 and T FO , the freezeout hypersurface changes significantly with c 2 s , with smaller c 2 s having freezeout hypersurface that spreads out further away from the origin (the propagation of the information that the system is expanding takes longer to propagate).Measurements of the speed of sound at midrapidity are thus less accurate for for a fixed T 0 and T FO , as it can be also remarked in fig. 1.The spatial spread-out of the hypersurface is confirmed in fig.3b, where we see a significant increase in the typical distance of the freezeout hypersurface when decreasing c 2 s from 1/3 to 1/4.Now we assess the effect of the width of initial temperature profile on the slope observable.To that end, we employ eq. ( 27) with σ = 5, 7.5, 10 fm, which are displayed in fig.4a.Similarly to the results of fig.1a, we employ fits to estimate d ln⟨p T ⟩/d ln N mr using globally-rescaled initial temperature profile, so that T 0 = αT ⋆ 0 , with α = 1.000, 1.005, • • • , 1.100 and T ⋆ 0 = 0.500 GeV.Figures 4b and 4c, portray results for c 2 s = 1/4 and 1/3, respectively.We see once again that the slope estimates are in general more accurate the smaller the value of T FO /T 0 is.Moreover, in both panels we see that the agreement is also improved the smaller the value of σ is, i.e., the more sharply peaked the initial temperature profile is.
As discussed in the previous sections, the slope estimate captures the value of c 2 s better if we take low values of the freezeout temperature T FO .In this regime, for Gubser flow, the constant temperature curves in the τ -r plane become highly peaked around the line τ = r, and the form of the freezeout hypersurface can be approximated by a simple function.For sufficiently low temperatures, the multiplicity and mean transverse momentum integrands in eq. ( 25) and eq.(26a) are dominated by the contributions stemming from the τ = r region, which we show in Appendix B to be given by7 3  26) with initial profiles given by eq. ( 30)) and the low T FO expansion formula (eq.( 33)).
We discuss in Appendix B why finite T FO corrections for N mr are larger than for N mr ⟨p T ⟩, hence the results are truncated at different orders in T FO /T 0 .
From the above expressions, we derive where the α derivatives are performed by considering the family of rescaled initial temperature profiles T 0 = αT ⋆ 0 discussed in Section II A. We see that the slope prescription indeed recovers the value of the speed of sound c 2 s = 1/3 for Gubser flow when T FO → 0. The result does not depend on the normalization constants entering in eq.(31).
The leading T FO /T 0 correction depends on the Gubser inverse width q and initial Milne time τ 0 through the factor (1 + q 2 τ 2 0 ).This shows clearly that the midrapidity multiplicity and average transverse momentum do retain a memory of the initial conditions through q and τ 0 , a memory which fades as T FO → 0. This dependence on the initial conditions is a reminder that treating the system as a thermodynamic one is only possible at sufficiently late times; eq. ( 32) provides a direct quantification of this effect.
Note that in general, the values usually employed for q and τ 0 in heavy-ion applications are such that qτ 0 is small.Equation (32) can thus be approximated well by the simpler expression In fig. 5, we assess the agreement between the full ultra-relativistic Cooper-Frye formulas (eq.( 26) with profiles given by eq. ( 30)) and expression (33) in the range 0 ≤ T FO /T 0 ≤ 0.2 GeV for two different values of the Gubser inverse width q.We see that eq. ( 33) captures very well the T FO /T 0 of the slope observable, with the numerical calculations (points) showing only slight deviations from the analytical formula (dashed lines).

V. SUMMARY
In this work, we studied the relation between the speed of sound of a fluid and its final particle multiplicity and average energy.Focusing on inviscid hydrodynamics with a constant speed of sound, we studied initial conditions where the temperature is scaled by an overall constant without change in its shape, a simplified model of ultracentral heavy-ion collisions.
Within this model, the average particle energy depends on the particle multiplicity strictly as a power law, with the power set by the square of the speed of sound, E tot /N tot ∝ N c 2 s tot , but only if the total energy and multiplicity of all particles can be measured.When restricting the set of measured final particles, the relation between the multiplicity, the energy and the speed of sound is only recovered for an infinitely long-lived fluid, which we characterized by an infinitely low freezeout temperature T FO → 0. Using three-dimensional hydrodynamic simulations, we confirmed the effect of rapidity cuts or transverse momentum cuts on the relation between the average energy E/N and the multiplicity N .We found a non-trivial dependence on rapidity due to the early freezeout of large-rapidity fluid elements when a fixed-temperature freezeout criteria is used.
To study in greater detail the midrapidity case, we used a simpler model with longitudinal boost invariance and cylindrical symmetry in the transverse plane.We focused on massless bosons produced at midrapidity, where the average energy corresponds to the average transverse momentum.We saw in this model that d ln⟨p T ⟩/d ln N mr converges as expected to c 2 s for asymptotically low freezeout temperatures.For a fixed ratio of the freezeout temperature to plasma's initial maximum temperature, T FO /T 0 , we saw that the convergence of d ln⟨p T ⟩/d ln N mr to c 2 s is faster for smaller c 2 s .Building on the above model, we used Gubser hydrodynamic solution to obtain expressions for d ln⟨p T ⟩/d ln N mr .The resulting formula quantifies explicitly the dependence of the latter quantity on the initial conditions and the freezeout temperature of the system, clarifying how and when the system approaches a "thermodynamic" limit where only conserved quantities are relevant.
While our model of ultra-relativistic collisions is a simplified one, our results suggest that the speed of sound will be best accessed with the highest possible center-of-mass energy collisions, where T FO /T 0 corrections will be smallest.Conversely, using eq.(1) in lower energy collisions may lead to inaccurate extractions of the speed of sound if the plasma's initial temperature is too low.Moreover, the non-trivial dependence on rapidity cuts observed in our model indicates future studies should be vigilant for this effect.
An important next step is to study temperature-dependent speeds of sound, which is the relevant scenario for nuclear matter.The effective temperature entering in eq. ( 1) is an essential ingredient of the original formulation of the slope observable in ref. [26], and the constant speed of sound model employed in the present work does not provide insights into the role of this parameter.In a different direction, it will be relevant to extend our reasoning to dissipative fluids and more general initial conditions.These extensions will depend more on numerical insights than analytical ones, and will be the subject of separate publications.In this Appendix8 , we shall provide further details on the derivation of eqs.( 25) and (26).In what follows, we parametrize the particle momentum p µ so that the transverse momentum 3-vector, ⃗ p T , points in the x-direction, thus (p µ ) t,x,ys,z = (m T cosh y, p T , 0, m T sinh y), where m T = p 2 T + m 2 , which reduces to m T = p T for massless particles, y is the particle momentumrapidity.In this case, the hypersurface element in cylindrical coordinates reads (dΣ µ ) t,r,ϕ,z = (rdrdϕdz, rdϕdzdt, 0, rdϕdrdt).For a cylindrically-symmetric, longitudinally boost-invariant flow, the fluid 4-vector can be parametrized as (u µ ) t,r,ϕ,z = (cosh κ cosh η s , sinh κ cosh η s , 0, sinh η s ).Then, we have that u µ p µ = m T cosh κ cosh(η s − y) − p T sinh κ cos ϕ.Besides, but, since we shall integrate η s ∈ (−∞, +∞), the second term of this expression does not contribute.Then, we have for the ultrarelativistic (m T ≈ p T ) particle yield at midrapidity (y = 0), where, from the first to the second line we have already performed the angular integral in the transverse momentum 3-vector, ⃗ p T , from the second to third line, we used the fact that From the third to the fourth line, we employ the variable transformation u = cos ϕ and perform the corresponding integrals.From the fourth to the fifth line, we perform the variable transformation v = sech 2 η s , which is performed in the step from the fifth to the sixth line.With analogous steps, we can also compute the mean transverse momentum as A2) where, from the second to third line, we used the fact that Appendix B: Gubser flow particlization and the low TFO limit For small values of the parameter q, the Gubser solution (30) for inviscid fluids exhibits initial temperature and flow profiles that share many similarities to the initial conditions expected in ultrarelativistic heavy-ion collisions: monotonically decreasing temperature profile in the radial direction with a peak at r = 0, and larger flow velocities in regions of larger initial temperature gradients.These features are illustrated in figs.6a and 6b for three different small values of q, at an initial time τ = τ 0 = 0.4 fm.
In fig.6c, we see, similarly with Section IV, that the low T FO hypersurfaces are sharply peaked.In the present case, this occurs exactly at τ = r.
In the T FO → 0 limit, the multiplicity N mr and transverse momentum N mr ⟨p T ⟩ are dominated by contributions originating from near the r = τ line.This can be seen when plotting the integrands of eqs.( 25) and (26a) in terms of the dimensionless coordinates X = q(τ − r), Y = q(τ + r), instead of τ and r, such that X near 0 corresponds to τ ≈ r.In fig. 7, we plot cosh κ (dr/dX)−sinh κ (dτ /dX) (the N mr integrand) and F (κ) cosh κ (dr/dX) − G(κ) sinh κ (dτ /dX) (the N mr ⟨p T ⟩ integrand) as a function of X, where κ = κ(X, Y (X)) is the fluid rapidity at the freezeout surface defined by Y (X), which is defined by implicitly the freezeout contour also shows that the N mr ⟨p T ⟩ integrand is much more sharply peaked around X = 0 than the N mr integrand.This will have an effect on how we perform the T FO → 0 expansion.In (X, Y ) coordinates, the multiplicity and mean transverse momentum integrals become The integration limits X ± are obtained from the original τ, r contour limits as follows: X + corresponds to the boundary at r = 0 and X − corresponds to the boundary at τ = τ 0 .More explicitly, X + is the real (positive) solution of the polynomial equation −4q 3 T 3 0 + T 3 FO X 5 + + 2T 3 FO X 3 + + T 3 FO X + = 0 and X − is the real (negative) solution of the polynomial equation τ In the limit where T FO → 0, we have that X − ≈ − √ 2[( T 3 0 q 2 )/(τ 0 T 3 FO )] 1/4 → −∞, X + ≈ 2 2/5 [( T 3 0 q 3 )/(T 3 FO )] 1/5 → +∞, Y (X) ≈ 2q T0 T FO 1 (1 + X 2 ) 1/3 → +∞, (B2) where the arrows denote the limit of the corresponding quantity in the strict limit as T FO → 0. Thus, this limit amounts to taking the integration boundaries to ±∞ and expanding the integrand in an asymptotic series as Y (X) → ∞.At leading order, the multiplicity integral becomes   (B4) Calculating O (T FO ) corrections to the multiplicity and transverse momentum is considerably more challenging.A key simplification in the eqs.(B3) and (B4) is the fact that the integrand is dominated by the X = 0 regime where Y is very large.Corrections of order O (T FO ) originate from regions where X is not small and Y is not large, making series expansion more challenging.On the other hand, as argued above with fig. 7, since the integrand of N mr ⟨p T ⟩ is more sharply peaked around X = 0, we do expect smaller O (T FO ) corrections; we verified this conclusion numerically.It is thus sufficient to find the O (T FO ) to the multiplicity, eq.(B3).
Guidance on the O (T FO ) to the multiplicity can be obtained from investigating corrections stemming from the fact that the integration interval is finite.One obtains The above expression is not a systematic expansion9 in T FO , but it does capture some of the dependence of the multiplicity on T FO .Based on this insight and on numerical tests, we found that the following expression captures the T FO dependence of N mr to very high accuracy: Agreement with numerical results is shown in fig.8, where eq.(B6) is labeled as "proposed NLO" and eq.(B5) is labeled as "partial NLO correction", both normalized with the LO expression (B3).In fact, in panel 8b, we see that the error of proposal (B6) is at most 0.2% in the interval considered.
FIG. 1: (a) Normalized average transverse momentum ⟨p T ⟩ as a function of normalized multiplicity N (see text for normalization prescription), both at midrapidity, for three constant speeds of sound c s , compared with the slope extracted from a linear fit of ln⟨p T ⟩ vs ln N ; (b) Dependence of the slope on the freezeout temperature T FO .

p∆η s / 2 TFIG. 2 :
FIG. 2: Slope observable as a function of size of rapidity window for different freezeout temperatures and p T cuts for c 2 s = 1/3 (a) and c 2 s = 1/4 (c).Slope estimation from energy and entropy flux through different isothermal hypersurfaces as a function of spatial rapidity window for c 2 s = 1/3 (b) and c 2 s = 1/4 (d).
resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231 using NERSC award NP-ERCAP0029957.Appendix A: Details for the particlization formulas