Smallest QCD Droplet and Multiparticle Correlations in pp Collisions

TL;DR

This paper investigates the applicability of hydrodynamics in proton-proton collisions, using the Gubser solution to model initial state fluctuations and predict multiparticle correlations, aligning well with experimental data and suggesting future tests at lower multiplicities.

Contribution

It introduces a model-independent approach combining initial state fluctuations with the Gubser solution to accurately describe multiparticle correlations in pp collisions.

Findings

Hydrodynamics is valid in certain pp collision regimes.

The model reproduces the sign and magnitude of multiparticle correlations.

Predicted sign change in correlations at lower multiplicities.

Abstract

The collective evolution of produced matter in heavy-ion collisions is effectively described by hydrodynamics from time scales greater than the inverse of the temperature, $\tau \gtrsim 1/T$. In the context of the Gubser solution, I show that the hydrodynamization condition $\tau \, T \gtrsim 1$ is translated into an allowed domain in the spatial system size and the final multiplicity for hydrodynamics applicability. It turns out that the flow measurements in pp collisions are inside the domain of validity. I predict that by approaching the boundaries of the allowed domain the hydrodynamic response to the initial ellipticity changes its sign. I follow a rather model-independent approach for the initial state where, instead of modeling the initial energy density of individual events, the initial system size and ellipticity event-by-event fluctuation are modeled. The model, initial state fluctuation+Gubser solution+Cooper-Frye freeze-out, describes the multiplicity and transverse momentum dependence of two-point and four-point correlation functions ($c_2\{2\}$ and $c_2\{4\}$) in an accurate agreement with pp collision experimental measurements. In particular, the sign of the four-point correlation function is the same as the observation, which failed to be described correctly in previous studies. I also predict a signal for the sign change in the hydrodynamic response that can be inspected in future experimental measurements of two-point and four-point correlation functions at lower multiplicities.

Tables

Table 1. Table 1: The rms radius fluctuation width and the kurtosis of ellipticity fluctuation from AMPT , and Γ 2 v superscript subscript Γ 2 𝑣 \Gamma_{2}^{v} from AMPT + GubsHyd .

${\sigma }_g$ [fm]	${\sigma }_r^{\text{AMPT}}$ [fm]	${\Gamma }_2^{ϵ\left(\text{AMPT}\right)}$	${\Gamma }_2^{v\left(\text{AMPT+Gubs}\right)}$
0.5	0.48	0.53	0.80
0.4	0.41	0.18	0.53
0.3	0.35	-0.17	0.26
0.2	0.30	-0.48	0.01
0.1	0.26	-0.73	-0.20

Equations

$${ϵ}_b(\rho ,\theta )=\frac{1}{{\tau }^4}\frac{{\widehat{ϵ}}_0}{{\cosh }^{8/3}\rho },$$

$${ϵ}_{\text{tot}}=\int rdrd\varphi ϵ(\tau ,r,\varphi ),$$

$$r_{\text{rms}}^2=\frac{1}{2{ϵ}_{\text{tot}}}\int rdrd\varphi r^2ϵ(\tau ,r,\varphi ),$$

$${ϵ}_2=-\frac{1}{2r_{\text{rms}}^2{ϵ}_{\text{tot}}}\int rdrd\varphi r^2\cos (2\varphi )ϵ(\tau ,r,\varphi ).$$

$${\widehat{ϵ}}_0=\frac{3{ϵ}_{\text{tot}}{r}_{\text{rms}}^2}{4\pi {\cosh }^{4/3}{\rho }_{\text{hyd}}}+\mathcal{O}({ϵ}_2),$$

$$\frac{1}{q^2}=r_{\text{rms}}^2\left(1+3{\tanh }^2{\rho }_{\text{hyd}}\right)+\mathcal{O}({ϵ}_2),$$

$$\lambda =\left(\sqrt{5\pi }/3\right){ϵ}_2+\mathcal{O}({ϵ}_2^2).$$

$${\rho }_{\text{hyd}}=-\mathrm{arccosh}\left[{\left(\frac{3r_{\text{rms}}^2{ϵ}_{\text{tot}}}{4\pi {\tau }_{\text{hyd}}^4{ϵ}_{\text{fo}}}\right)}^{1/4}\right]+\mathcal{O}({ϵ}_2).$$

$$\frac{dN}{d\mathbf{p}}=-\frac{g}{(2\pi {)}^3}\int p^{\mu }d{\Sigma }_{\mu }\exp \left[p_{\mu }{u}^{\mu }/T_{\text{fo}}\right],$$

$$\begin{array}{rl}& \cos {\theta }_{\text{fo}}(\rho )=\tanh \rho +\frac{1}{q}{\left(\frac{{ϵ}_{\text{fo}}}{\widehat{{ϵ}_0}}\right)}^{1/4}{\mathrm{sech}}^{1/3}\rho .\end{array}$$

$$r_{\text{rms}}\ge r_{\text{crit}},r_{\text{crit}}=\sqrt{\frac{4\pi }{3}}\left[{\tau }_{\text{hyd}}^2{ϵ}_{\text{fo}}^{1/2}\right]{ϵ}_{\text{tot}}^{-1/2}.$$

$$r_{\text{rms}}\le R_{\text{crit}}({ϵ}_{\text{tot}}),$$

$$v_2(n_{\text{tot}},r_{\text{rms}})=k_2(n_{\text{tot}},r_{\text{rms}}){ϵ}_2.$$

$$v_2=\chi k_2({ϵ}_{\text{tot}},r_{\text{rms}}){ϵ}_2.$$

$$c_2\{2\}$$

$$c_2\{4\}$$

$$v_2^{\text{Gub}}\{2k\}=\chi v_2^{\text{IdealGub}}\{2k\}.$$

$$v_2\approx {\tilde{k}}_2{ϵ}_2,$$

$$v_2^{\text{Glaub}}\{2k\}={\tilde{k}}_2{ϵ}_2\{2k\}.$$

$$p_{\text{init}}({ϵ}_2,r_{\text{rms}})=p_{ϵ}({ϵ}_2)p_r(r_{\text{rms}}).$$

$$k_2^{\text{caus}}\equiv \chi k_2(n_{\text{tot}},r_{\text{rms}}).$$

$$\begin{array}{r}{p}_v(v_2;n_{\text{tot}})=\int \frac{d{r}_{\text{rms}}}{k_2^{\text{caus}}}{p}_{ϵ}(v_2/k_2^{\text{caus}})p_r(r_{\text{rms}}),\end{array}$$

$$p_{ϵ}({ϵ}_2)=\frac{{ϵ}_2}{{\sigma }_{ϵ}^2}{e}^{-{ϵ}_2^2/2{\sigma }_{ϵ}^2}[1+\frac{{\Gamma }_2^{ϵ}}{2}{L}_2({ϵ}_2^2/2{\sigma }_{ϵ}^2)+\cdots ],$$

$${\Gamma }_2^{ϵ}\equiv -({ϵ}_2\{4\}/{ϵ}_2\{2\}{)}^4$$

$$p_r(r_{\text{rms}})=\frac{r_{\text{rms}}}{{\sigma }_r^2}{e}^{-r_{\text{rms}}^2/2{\sigma }_r^2},$$

$$v_2\{2\}$$

$$v_2\{4\}$$

$$\chi {\sigma }_{ϵ}\approx 0.097,{\Gamma }_2^{ϵ}\approx -0.75,{\sigma }_r\approx 0.4[\text{fm}],$$

$$⟨N_{\text{ch}}⟩$$

$$\times N_{\text{tot}}(n_{\text{tot}},r_{\text{rms}},p_{T,\text{min}},p_{T,\text{max}},{\eta }_{\text{min}},{\eta }_{\text{max}})$$

$${\Gamma }_2^v=c_2\{4\}/c_2^2\{2\},$$

$$-0.39\lesssim {\Gamma }_2^{v(\text{ATLAS})}\lesssim -\mathrm{0.31.}$$

$${\Gamma }_2^{ϵ(\text{TREN})}\approx 0.09,{\sigma }_r^{\text{TREN}}\approx 0.47\text{[fm]}.$$

$$\widehat{ϵ}=\frac{{\widehat{ϵ}}_0}{{\cosh }^{8/3}\rho }\left[1+4\lambda {\delta }_2(\rho )y(\theta ,\varphi )\right]+\mathcal{O}({\lambda }^2),$$

Full text

I Introduction

In 2010, CMS collaboration revealed a peculiar observation of long-range correlation in p–p collisions Khachatryan:2010gv which is considered as a signature of collective evolution. Later on, this observation was confirmed by different experimental collaborations for different small systems (p–p, pAu, dAu, 3HeAu, pPb) at LHC Abelev:2014mda ; Khachatryan:2016txc ; Aaboud:2016yar ; Aaboud:2017blb and RHIC Adamczyk:2015xjc ; PHENIX:2018lia . Over the past years, there have been ongoing debates on the origin of the observed correlation. Efforts to explain the observed phenomena have been made from different perspectives, e.g., to describe the phenomenon via kinetic theory Kurkela:2018qeb ; Kurkela:2019kip or to link the correlation to the initial stages of the collision (for review see Ref. Strickland:2018exs ).

The present paper belongs to the category of studies that intend to demonstrate the observed phenomena using conventional hydrodynamics Bozek:2011if ; Bozek:2013uha ; Niemi:2014wta ; Weller:2017tsr ; Mantysaari:2017cni ; Gallmeister:2018mcn ; Heinz:2019dbd . The strategy that I pursue in the present paper is the following: studying hydrodynamic evolution such that, first, it has essential features to explain the real data, and, second, it is still simple enough to monitor an event evolution anatomy clearly. To this end, the best choice in the author’s opinion is the analytical solution of relativistic hydrodynamic equations for conformal fluids, the Gubser solution, and perturbation on top of that Gubser:2010ze ; Gubser:2010ui . Despite the idealized assumptions for the symmetries of the system, the Gubser solution has been used to investigate, to an extent, realistic scenarios in heavy-ion collisions. As an example, the power spectrum of heavy-ion collisions is obtained by studying the evolution of narrow peaked hot spots on top of a smooth background in Refs. Staig:2011wj ; Gorda:2014msa .

In the present paper, I consider the rotationally symmetric Gubser solution, which is perturbed to acquire an initial elliptical shape Gubser:2010ui . The free parameters of the Gubser solution are carefully mapped to the more standard quantities such as the root mean square (rms) radius, ellipticity, and the total energy in the transverse direction of the system. Using this map together with the equation of state (e.o.s.), the initiation time of the evolution, the mass of the final particles, and the freeze-out temperature, one obtains an estimation for the elliptic flow of any given initial energy distribution. Considering the consistency of the Gubser solution, I find the domains of hydrodynamics applicability in terms of the spatial system size and final multiplicity. One observes that flow measurements of p–p collisions are inside the domain of validity. I also describe the measured two-particle, $c_{2}\{2\}$, and four-particle, $c_{2}\{4\}$, correlations as a function of transverse momentum and charge multiplicity. Moreover, I predict an experimental signal in multiparticle correlation functions in p–p collisions to indicate whether one approaches the boundaries of the validity domain. Indeed, due to the several simplifications and idealizations in the Gubser solution, the present model is not applicable to gain insight into hydrodynamic transport coefficients such as shear viscosity over entropy density. A more realistic hydrodynamic simulation is essential for such studies. I justify, however, that the computations are accurate enough for the purposes I pursue in this paper.

The paper is organized as follows: In section II, I investigate the domain of applicability of hydrodynamics as a function of the initial system size and the total transverse energy. In section III, the model is introduced and compared with the outcome with iEBE-VISHNU. The model is compared with experimental p–p collision measurements in section IV, and finally one finds the conclusion in section V. The computational details of the model are presented in appendix A. Further comparisons between predictions of the model and iEBE-VISHNU can be found in appendix B.

II Smallest system Size for Hydrodynamic Description

I start by highlighting the main aspects of the Gubser solution Gubser:2010ze ; Gubser:2010ui . Before that, I determine the coordinate systems employed in this paper. I interchangeably use hyperbolic-cylindrical coordinates $x^{\mu}=(\tau,r,\phi,\eta)$ ($\tau=\sqrt{t^{2}-z^{2}}$ is proper time, $\eta=\operatorname{arctanh}(z/t)$ is space-time pseudorapidity, and $(r,\phi)$ is the polar coordinate in the transverse plane) and the de Sitter coordinate $\hat{x}^{\mu}=(\rho,\theta,\phi,\eta)$ ($\rho$ is de Sitter time, and $\theta$ is de Sitter radial coordinate in the transverse direction). Two coordinates are related to each other via $q\,\tau=\operatorname{sech}\rho/(\cos\theta-\tanh\rho)$ and $q\,r=\sin\theta/(\cos\theta-\tanh\rho)$ where $q$ is a free parameter.

In the de Sitter coordinate, the Gubser solution for the energy density of a boost-invariant ideal conformal fluid reads as

[TABLE]

where $\tau$ should be written in terms of $\theta$ and $\rho$ with the transformation mentioned earlier. In the above, $\hat{\epsilon}_{0}$ is a free parameter. An elliptic perturbation on top of the rotationally symmetric background solution (i.e. Eq. (1)) leads to $\epsilon(\rho,\theta,\phi)\approx\epsilon_{b}(\rho,\theta)\left[1+4\lambda\,\delta_{2}(\rho)\,y(\theta,\phi)\right]$ with $y(\theta,\phi)=-\sqrt{3/8}Y_{2,2}+Y_{2,0}/2-\sqrt{3/8}Y_{2,-2}$ ($Y_{l,m}\equiv Y_{l,m}(\theta,\phi)$ is the spherical harmonic function). The fluid velocity is given by $u_{\mu}=\tau(\partial\hat{x}^{\nu}/\partial x^{\mu})\hat{u}_{\nu}$ where $\hat{u}_{\mu}=(-1,\;\lambda\,\delta\hat{u}_{\theta},\;\lambda\,\delta\hat{u}_{\phi},0)$, $\lambda$ is a small free parameter, and $\delta\hat{u}_{i}=\tilde{\delta}_{2}(\rho)\;\partial_{i}y(\theta,\phi)$ for $i=\theta,\phi$ Gubser:2010ze ; Gubser:2010ui ; Staig:2011wj (see a more detailed review in appendix A). The hydrodynamic equations lead to a linear system of differential equations for $\delta_{2}(\rho)$ and $\tilde{\delta}_{2}(\rho)$ which can be solved analytically for ideal hydrodynamics. The isotropic initial fluid velocity assumption at de Sitter time $\rho=\rho_{\text{hyd}}$ fixes the initial value for the system of equations as $\delta_{2}(\rho_{\text{hyd}})=1$ and $\tilde{\delta}_{2}(\rho_{\text{hyd}})=0$. In this paper, I do not present the Navier-Stokes solution, which reveals instability as one expects. The problem could be cured in a causal hydrodynamics framework. However, the Gubser solution for causal hydrodynamics has been found only in certain limits Marrochio:2013wla ; Pang:2014ipa , and perturbation on top of causal Gubser flow has not been done yet. In addition to the ignored dissipative effects, the conformal symmetry prevents one from choosing a realistic e.o.s. Moreover, the hydrodynamic solution is limited to an elliptic perturbation (a linear hydrodynamic response) which also means that interaction between the ellipticity mode with itself (e.g., cubic response Noronha-Hostler:2015dbi ) and other mode-mode interactions are ignored. I address the effect of all these idealizations before I start model/data comparison.

From now on, as an abbreviation, the Gubser hydrodynamic solution together with Cooper-Frye prescription Cooper:1974mv for the freeze-out is called GubsHyd. A detail of the “model” can be found in appendix A. A difference between GubsHyd and the conventional numerical hydrodynamic solutions is in their initiations. In the numerical hydrodynamic computations, the initial energy density is prepared on a $\tau=\text{cte}$ surface (Fig. 1 (left)) while GubsHyd is initiated on a surface with condition $\rho=\rho_{\text{hyd}}=\text{cte}$ (Fig. 1 (right)). Referring to Eq. (1), one simply finds that $\rho=\text{cte}$ is equivalent to $\tau\epsilon^{1/4}=\text{cte}$ which is subsequently equivalent to $\tau T=\text{cte}$ by using e.o.s. $\epsilon=C_{0}T^{4}$.

It is worth reminding the reader that although I explicitly present the ideal Gubser solution in Eq. (1), I implicitly assume it is only an approximation of a more general causal hydrodynamics in which the system is not in perfect local equilibrium. As a result, to apply hydrodynamics to a far-from-equilibrium initial state, a finite time is needed for hydrodynamization. In this context, the $\tau T=\text{cte}$ surface can be understood as follows: regarding studies about hydrodynamization, specifically, the computations from gauge/gravity duality, the evolution of a boost-invariant system with translational and rotational symmetry in the transverse space (Bjorken symmetry) is attracted to the hydrodynamic solutions after the time $\tau\sim 1/T$ Chesler:2009cy ; Heller:2011ju . It has been demonstrated that the hydrodynamic gradient expansion ($1/\tau T$ expansion in this context) is divergent, and its divergence is due to the presence of nonanalytic contributions $\propto e^{-z_{0}\tau T}$ known as nonhydrodynamic modes Heller:2013fn ; Heller:2015dha ; Romatschke:2017vte ; Romatschke:2017acs ; Heller:2016rtz . Here, the coefficient $z_{0}\sim\mathcal{O}(1)$ is a positive real numerical factor depending on the underlying microscopic theory. Therefore, the nonhydrodynamic modes die out at $\tau T\sim 1$, and the hydrodynamic description works appropriately afterward.111The attractor solution has been observed in an expanding ultrarelativistic gas of hard spheres with Bjorken symmetry. It has been shown that the causal hydrodynamics can describe this system even when gradients are significant at all times of evolution Denicol:2019lio .

The attractor solution of a system that goes through Gubser flow has been studied in Ref. Behtash:2017wqg within the relativistic kinetic theory framework (see also Ref. Chattopadhyay:2018apf ). This paper shows that hydrodynamic quantities which are parametrized with variable $w=\tanh\rho/\tau T$ approach to an attractor at a specific value of $w$, namely $w\sim w_{0}$. Since $\tau T\sim\tau\epsilon^{1/4}\sim\cosh^{-2/3}\rho$, the condition $w\sim w_{0}$ is equivalent to the condition $\tau T\sim\text{cte}$, similar to what has been found for the systems with Bjorken symmetry. In Gubser flow, however, the temperature drops when one moves from the center to the tail of the initial energy density. As a result, it is plausible to consider that the hydrodynamization happens at different proper times depending on the transverse radius. This argument allows one to define the hydrodynamization surface, the surface on which $\tau\,T(\tau,r)$ is constant and is in the order of unity [see Fig. 1 (right)]. For a system with fewer symmetries and more complicated initial conditions, one would expect that the nonhydrodynamic modes decay at $\tau\sim 1/T(\tau,\vec{x})$, leading to a more complicated hydrodynamization surface $\tau=\tau(\vec{x})$.

To compare the GubsHyd predictions with real experimental data, one translates the free parameters of the Gubser solution $(q,\hat{\epsilon}_{0},\lambda,\rho_{\text{hyd}})$ to more standard quantities: total transverse energy $\epsilon_{\text{tot}}$, rms radius $r_{\text{rms}}$, and ellipticity $\epsilon_{2}$,

[TABLE]

The measure of these integrals reads as $\tau^{2}\cosh^{2}\rho_{\text{hyd}}\,\sin\theta\,d\theta d\phi$ in the de Sitter coordinates. Using the measure together with the fact that the solution is initiated on the $\rho=\rho_{\text{hyd}}$ surface, one obtains

[TABLE]

For finding $\rho_{\text{hyd}}$ in terms of standard quantities, one notes that some parts of the initial energy density are about to freeze-out immediately after the initiation because the energy density of these parts is equal to the freeze-out energy density $\epsilon_{\text{fo}}$. As a result, at the hydrodynamization surface, one has $\tau\epsilon^{1/4}=\tau_{\text{hyd}}\epsilon_{\text{fo}}^{1/4}$. Subsequently, by employing Eq. (1) and substituting $\hat{\epsilon}_{0}$ in terms of standard quantities (Eq. (3a)), one finds

[TABLE]

I employ Cooper-Frye prescription Cooper:1974mv with Boltzmann equilibrium distribution to obtain the associated final particle distribution function:

[TABLE]

where $d\boldsymbol{p}\equiv d^{3}p/E$ and $\Sigma^{\mu}=(\rho,\theta_{\text{fo}}(\rho,\phi),\phi,\eta)$ indicates the freeze-out surface specified by equation $\epsilon(\rho,\theta_{\text{fo}},\phi)=\epsilon_{\text{fo}}$. In Eq. (5), $g$ is the degeneracy constant and $T_{\text{fo}}$ is the freeze-out temperature related to $\epsilon_{\text{fo}}$ via e.o.s. The equation $\epsilon(\rho,\theta_{\text{fo}},\phi)=\epsilon_{\text{fo}}$ for the unperturbed Gubser solution (Eq. (1)) leads to the following equation for the freeze-out surface:

[TABLE]

Red curves show the freeze-out surfaces in Fig. 2 for systems with three different sizes.

Initiating the evolution on the hydrodynamization surface leads to an interesting conclusion. One defines $\bar{\tau}_{\text{fo}}$ as the time when the last fluid cell of the system is frozen out (blue open circles in Fig. 2). On the other hand, by definition, $\tau_{\text{hyd}}$ on the hydrodynamization surface is also a member of the freeze-out surface (blue bullets show $\tau_{\text{hyd}}$ location in Fig. 2). One can interpret $\tau_{\text{hyd}}$ as the time at which the last fluid cell is hydrodynamized. One notes that the condition $\tau_{\text{hyd}}\leq\bar{\tau}_{\text{fo}}$ should always be satisfied; otherwise, the last frozen-out cell has not enough time to be hydrodynamized. It turns out that there are some values for $(r_{\text{rms}},\epsilon_{\text{tot}},\epsilon_{\text{fo}},\tau_{\text{hyd}})$ that correspond to no valid hydrodynamic solution for making a hydrodynamization surface. With that, one concludes a lower bound on the system size. Considering that the argument in Eq. (4) should be greater than unity to have a real-valued $\rho_{\text{hyd}}$, one obtains the following lower bound on the system size,

[TABLE]

In Fig. 2, the freeze-out surface for a system at a critical size is depicted in the right panel.

The above arguments can be considered as a hint for the presence of the same lower bound for any nonideal/nonconformal hydrodynamic systems. Naïve reasoning is as follows: Consider a nonequilibrium system with a fixed amount of total transverse energy $\epsilon_{\text{tot}}$. The smaller the system size, the sharper the concentration of energy. Compared to a system with a larger size, a system with more concentration of energy would expand faster because of the more considerable “pressure” gradient. Therefore, one might consider cases in which the system expands so fast that the nonhydrodynamic modes do not have enough time to decay before the nonequilibrium system reaches the phase transition boundary. For such a system, equation $\tau\epsilon^{1/4}(\tau,\vec{x})=\tau_{\text{hyd}}\epsilon_{\text{fo}}^{1/4}$ has no solution, and the hydrodynamization surface $\tau(\vec{x})$ does not exist. For more robust evidence, computations beyond the hydrodynamic regime are needed. In fact, the bound (7) is compatible with that obtained from studying two colliding shock waves in the context of numerical holography Chesler:2015bba ; Chesler:2016ceu . One can estimate the averaged initial energy density as $\bar{\epsilon}_{\text{init}}\sim\epsilon_{\text{tot}}/\pi r_{\text{crit}}^{2}$. Using $T_{\text{eff}}=(4\bar{\epsilon}_{\text{init}}/3\pi^{4})^{1/4}$ (Eq. (7) in Ref. Chesler:2015bba ), the condition $r_{\text{crit}}\,T_{\text{eff}}\sim 1$ (Eq. (10) in Ref. Chesler:2015bba ) is equivalent to $r_{\text{crit}}\,\epsilon_{\text{tot}}^{1/2}\sim 1$ in Eq. (7). The other possibility is using the kinetic theory with nontrivial energy distribution in the transverse direction, such as that studied in Refs. Kurkela:2018qeb ; Kurkela:2019kip . Studying the hydrodynamic attractors for such a system would lead to evidence for the presence of a lower bound on the system size at which hydrodynamics is applicable.

From Eq. (6), one can also find the de Sitter time when the last fluid cell is frozen out via $\cos\theta_{\text{fo}}(\rho_{\text{max}})=1$. Consequently, the condition $\rho_{\text{hyd}}=\rho_{\text{max}}$ means there is no “time” left for hydrodynamic evolution. With that, one concludes an upper bound,

[TABLE]

for the system size corresponding to the case that there is not enough energy density deposited into the given region to produce a deconfined matter. Translating variables into standard quantities, one can solve the equation $\cos\theta_{\text{fo}}(\rho_{\text{hyd}})=1$ in Eq. (6) to find $R_{\text{crit}}(\epsilon_{\text{tot}})$ in Eq. (8).

III Hydrodynamic response to initial ellipticity

The domain of hydrodynamic validity has been investigated in the previous section. To determine whether the experimental measurements for p–p collision are inside this domain, one needs to relate the hydrodynamic solutions to the final particle distribution. Having found the freeze-out surface, one can compute the final particle distribution by performing the integrals in the Cooper-Frye formula. The details of the computation are presented in appendix A. Furthermore, these computations are implemented in a Mathematica package available as an ancillary file. In the following, I briefly mention some of the functions discussed in appendix A: the total multiplicity in the unit rapidity $n_{\text{tot}}(r_{\text{rms}},\epsilon_{\text{tot}})$ ($n_{\text{tot}}\equiv dN/dy_{p}$); the translation of critical values $r_{\text{crit}}(\epsilon_{\text{tot}})$ and $R_{\text{crit}}(\epsilon_{\text{tot}})$ in terms of multiplicity in unit rapidity, $n_{\text{crit}}(r_{\text{rms}})$ and $N_{\text{crit}}(r_{\text{rms}})$; the particle distributions in unit space-time pseudorapidity $dN/p_{T}dp_{T}d\eta$ and $dN/d\eta$; and finally and most importantly $k_{2}(n_{\text{tot}},r_{\text{rms}})$ and $k_{2}(\epsilon_{\text{tot}},r_{\text{rms}})$ (and their $p_{T}$ dependent versions) which relate the initial ellipticity $\epsilon_{2}$ to the elliptic flow, the second Fourier coefficient of the final particle distribution in the azimuthal direction,

[TABLE]

The initial hydrodynamization time $\tau_{\text{hyd}}$, the final particle mass $m$, the freeze-out energy density $\epsilon_{\text{fo}}$, and the constant $C_{0}$ in e.o.s. are considered to be fixed in the model.

In Fig. 3, the contour plot of $k_{2}(n_{\text{tot}},r_{\text{rms}})$ for pions $m\approx 136\;$MeV is depicted for $0.3<p_{T}<3\,\text{GeV/c}$. Here, I have fixed $\tau_{\text{hyd}}=0.9\,\text{fm/c}$, $\epsilon_{\text{fo}}=0.18\;\text{GeV/fm}^{3}$, and $C_{0}=13$. The latter is chosen close to the plateau in s95p-v1 Huovinen:2009yb at temperatures above the critical point where most of the evolution period is in this range. The functions $n_{\text{crit}}(r_{\text{rms}})$ and $N_{\text{crit}}(r_{\text{rms}})$ are also shown with thick black curves. What one sees from the figure is that the system is hydrodynamized in the range $n_{\text{tot}}>n_{\text{crit}}(r_{\text{rms}})$. Moreover, in the valid domain of hydrodynamics, there are two different regions: the in-plane response region where particles are mostly emitted along the minor axis of the initial elliptic shape (blue contours) and the out-of-plane response regions where most of the particles are emitted along the major axis (red contours). As it will be discussed later, the signature of in-plane/out-of-plane transition is traceable in multiparticle correlations in p–p collisions.

Due to the idealizations that have been made, one might be doubtful about the applicability of the presented hydrodynamic model in explaining the experimental data. To show that GubsHyd is accurate enough for the purpose I pursue in this paper, I compare a causal hydrodynamic simulation, iEBE-VISHNU (MC-Glauber+VISH2+1) Shen:2014vra with GubsHyd. Since the interested quantities are those which are averaged over an ensemble of events, I show that all the corrections can be encapsulated in a generic constant $\chi$,

[TABLE]

It will be shown that $\chi$ can be absorbed into the ellipticity fluctuation width, a free parameter in the model, and has no impact on the predictions in the present paper.

I employ iEBE-VISHNU with MC-Glauber as an initial state model for Pb–Pb collisions at $\sqrt{s_{NN}}=2.76\,\text{TeV}$. The simulation is based on solving 2+1 boost-invariant causal hydrodynamic equations with a fixed $\eta/s=0.08$. The initial time of hydrodynamic evolution is $\tau=0.6\;\text{fm/c}$, and freeze-out energy density is fixed to $\epsilon_{\text{fo}}=0.18\;\text{GeV/fm}^{3}$. The afterburner is not included in this simulation. The events are classified into 16 bins between [math] and $80\%$ centrality classes, and at each bin 14 000 events are generated. The elliptic flow is computed for charged pions with the momentum in the range $0.28<p_{T}<4\,\text{GeV}$ at each event. After that, the probability density function (p.d.f.) of elliptic flow fluctuation $p_{v}(v_{2})$ at each centrality class is obtained, and eventually, I compute the two first cumulants of the p.d.f. Borghini:2001vi

[TABLE]

where $\langle\cdots\rangle_{v}$ stands for averaging with respect to $p_{v}(v_{2})$. Now, from the same set of initial states, I compute $(\epsilon_{\text{tot}},r_{\text{rms}},\epsilon_{2})$ at each event. Setting $\tau_{\text{hyd}}=0.6\,\text{fm/c}$, $m=136\;\text{MeV}$, $\epsilon_{\text{fo}}=0.18\;\text{GeV/fm}^{3}$, and $C_{0}=13$, I employ GubsHyd (the function in the second line of Eq. (LABEL:50)) for $0.28<p_{T}<4\,\text{GeV}$ and use $v_{2}=k_{2}(\epsilon_{\text{tot}},r_{\text{rms}})\,\epsilon_{2}$ to obtain $v_{2}$ at each event. Eventually, one has 14 000 events for each centrality class where its elliptic flow is evaluated from GubsHyd. Similar to iEBE-VISHNU events, one can compute $v_{2}^{\text{\tiny IdealGub}}\{2\}$ and $v_{2}^{\text{\tiny IdealGub}}\{4\}$ from GubsHyd as well. One notes from Eq.  (11) that the constant $\chi$ in Eq. (10) corrects the cumulants as

[TABLE]

The results are shown by red filled circles in Fig. 4 where $\chi=0.39$ is chosen. As seen from the figure, there is a good agreement between iEBE-VISHNU and GubsHyd in the given range of centralities.

Eq. (10) is different from that mentioned in Ref. Gardim:2011xv as a linear hydrodynamic response:

[TABLE]

where $\tilde{k}_{2}$ is a constant unlike $k_{2}(\epsilon_{\text{tot}},r_{\text{rms}})$. To show that the observed agreement is not trivially inherited from the initial state, one can obtain the same cumulants as Eq. (12) by employing Eq. (13):

[TABLE]

In the above, $\epsilon_{2}\{2k\}$ can be obtained from Eq. (11) by replacing the variable $v_{2}$ with $\epsilon_{2}$ while the average $\langle\cdots\rangle_{\epsilon}$ is performed using the distribution $p_{\epsilon}^{\text{\tiny Glaub}}(\epsilon)$ associated with the ellipticity fluctuation in the MC-Glauber model. In Fig. 4, $v_{2}^{\text{\tiny Glaub}}\{2\}$ and $v_{2}^{\text{\tiny Glaub}}\{4\}$ with $\tilde{k}_{2}=0.25$ are plotted by the green filled triangles. Although linear response explains the causal hydrodynamic simulations in central collisions, it fails to explain it in peripheral collisions.

The above investigation indicates that, concerning many event averages, the effect of idealizations in the GubsHyd model can be corrected approximately via a constant value $\chi$. The numerical value of $\chi$ is irrelevant in the analysis I perform in the next section. For completeness, the GubsHyd predictions have been compared with iEBE-VISHNU for $p_{T}$ spectrum, differential $v_{2}\{2\}$, and $v_{2}\{4\}$; and centrality dependence of the multiplicity for charged pions in appendix B.

IV Two- and Four-particle correlations in p–p collisions

The initial state is the final piece in the model presented in this paper. In this paper, I follow a rather model-independent approach for the initial state. Instead of generating an ensemble of initial energy densities, $\epsilon_{2}$ and $r_{\text{rms}}$ event-by-event fluctuations are modeled. By comparing the model with data, one finds the properties of the $\epsilon_{2}$ and $r_{\text{rms}}$ fluctuations. In the end, one can check different initial state models (or one model with different parameter tuning) to determine which one can reveal the same fluctuating properties.

Assuming $r_{\text{rms}}$ and $\epsilon_{2}$ fluctuate event by event independently, the p.d.f. of the two-dimensional (2D) variable $(\epsilon_{2},r_{\text{rms}})$ can be written as a product of ellipticity distribution $p_{\epsilon}(\epsilon_{2})$ and rms radius distribution $p_{r}(r_{\text{rms}})$;

[TABLE]

To find the p.d.f. of elliptic flow, one changes the variable $(\epsilon_{2},r_{\text{rms}})$ to $(v_{2},r_{\text{rms}})$ via Eq. (9). The new p.d.f. reads as $p_{\epsilon}(v_{2}/k_{2}^{\text{caus}})p_{r}(r_{\text{rms}})/k_{2}^{\text{caus}}$ where

[TABLE]

By averaging out the variable $r_{\text{rms}}$, a p.d.f. for $v_{2}$ fluctuation is obtained,

[TABLE]

where $n_{\text{tot}}$ in the argument appears as a parameter.

Now, I estimate each $p_{\epsilon}(\epsilon_{2})$ and $p_{r}(r_{\text{rms}})$ separately. The ellipticity distribution can be expanded as the following series Abbasi:2017ajp ; Mehrabpour:2018kjs :

[TABLE]

where $L_{2}(\epsilon_{2}^{2}/2\sigma_{\epsilon}^{2})$ is the second Laguerre polynomial and

[TABLE]

is the kurtosis of the distribution. One could use the power distribution introduced in Ref. Yan:2013laa for the ellipticity fluctuation in small systems. Compared to power distribution which is obtained for point-like sources, the distribution in Eq. (16) is more general. It turns out that the ellipsis in Eq. (16) expansion is irrelevant to $v_{2}\{2\}$ and $v_{2}\{4\}$ quantities. Therefore, concerning the ellipticity fluctuation, no physical assumption is imposed here. For $r_{\text{rms}}$, I assume a Gaussian distribution

[TABLE]

with free parameter $\sigma_{r}$ and ignore other corrections to it.

By substituting Eq. (16) into Eq. (15) and using relations in Eqs. (11) and $k_{2}^{\text{caus}}\approx\chi\,k_{2}$, one finally obtains the model predictions:

[TABLE]

with free parameters $\sigma_{\epsilon}$, $\Gamma_{2}^{\epsilon}$, and $\sigma_{r}$. In the above, $\langle\cdots\rangle_{r}$ refers to averaging with respect to $p_{r}(r_{\text{rms}})$. One obtains the $p_{T}$ dependent or $p_{T}$ integrated predictions for $v_{2}\{2\}$ and $v_{2}\{4\}$ from GubsHyd by substituting the $p_{T}$ dependent or $p_{T}$ integrated $k_{2}$ function into relations (19) (see Eq. (LABEL:a38)).

After adding the final piece into the model, one can compare it with experimental data. In Fig. 5, the measured $v_{2}\{2\}$ and $v_{2}\{4\}$ by CMS Khachatryan:2016txc and ATLAS Aaboud:2017blb collaborations for p–p collisions at $\sqrt{s}=13\;\text{TeV}$ are presented. The GubsHyd predictions for $v_{2}\{2\}$ (blue curve) and $v_{2}\{4\}$ (red dashed curve) are also shown where I have used the following parameters,

[TABLE]

in Eq. (19). To compare with experimental data, the parameter $n_{\text{tot}}$ is translated to the average charge multiplicity $\langle N_{\text{ch}}\rangle$ in the range $|\eta|<2.4$ via

[TABLE]

where $N_{\text{tot}}$ can be found in Eq. (LABEL:A22). In measuring $v_{2}\{4\}$ in the ATLAS collaboration, the three-subevent method is employed to reduce the non-flow effect Aaboud:2017blb . For that reason, I compare the model prediction for $v_{2}\{4\}$ only with ATLAS results. I did not present multiplicity dependence of $v_{2}\{2\}$ from the ATLAS collaboration (with peripheral subtraction method) because it follows the same trend as the CMS result (with $\eta$-gap $\Delta\eta>2$) and makes the plot hard to read.222The data points of ATLAS measurements for $v_{2}\{2\}$ are not available online. For the comparison with CMS data, I extracted the ATLAS measurement data points directly from the published figure. This is a second reason I have not used the ATLAS data in the figure. Given that the idealized model GubsHyd has been employed to extract the initial state fluctuation properties, the numbers in Eq. (20) should be assumed as an approximate estimation.

By comparing the model with experimental data, one could extract the parameters $\Gamma_{2}^{\epsilon}$ and $\sigma_{r}$. However, due to the lack of knowledge about constant $\chi$, the parameter $\sigma_{\epsilon}$ remains undetermined. The term $\chi\,\sigma_{\epsilon}$ is an overall factor in Eq. (19) and has no impact on multiplicity and $p_{T}$ dependence of $v_{2}\{2\}$ and $v_{2}\{4\}$. It means irrespective to the value of $\chi$ the agreement between GubsHyd prediction and experimental data for $v_{2}\{2\}$ and $v_{2}\{4\}$ is a triumph for the model. Moreover, an interesting physics prediction from the model can be tested by pushing the experimental measurements close to the critical multiplicity in the future. In Fig. 5 (left), the model predicts a valley around $\langle N_{\text{ch}}\rangle=10$ shown by a vertical dashed line. This is an indication of in-plane/out-of-plane transition. By translating $\langle N_{\text{ch}}\rangle$ to $n_{\text{tot}}$, a similar vertical dashed line is shown in Fig. 3 and a horizontal dashed line that locates $\langle r_{\text{rms}}\rangle_{r}$. The latter is shown as an indicator for the average value of the system size. Moving from right to left along the horizontal dashed line in Fig. 3, $k_{2}(n_{\text{tot}},r_{\text{rms}})$ changes sign from positive to negative. Since there are only even powers of $k_{2}$ in Eq. (19), the sign change makes a valley in $v_{2}\{2\}$ and $v_{2}\{4\}$ versus multiplicity. Indeed, nonflow effects are more pronounced at lower multiplicities, and observing the flow signals at lower multiplicity would be experimentally challenging. The discontinuity in $v_{2}\{4\}$ prediction is due to the fact that in a region of multiplicity, the terms in the bracket in Eq. (19b) turn to negative values. One notes that the predicted valley in multiplicity dependence of $v_{2}\{2\}$ and $v_{2}\{4\}$ is based on an idealized hydrodynamic model. A framework beyond hydrodynamics is needed for a more accurate prediction.

As a final remark, I comment on the sign of $c_{2}\{4\}=-(v_{2}\{4\})^{2}$ when one initiates the energy density with a realistic model. This model must contain nucleon substructure to produce enough fluctuation. To this end, the AMPT event generator Lin:2004en and T${}_{\text{R}}$ENTo with nucleon substructure Moreland:2014oya ; Moreland:2018gsh are examined as initial state models.

I follow Ref. Xu:2016hmp and employ the AMPT event generator in the string melting mode to generate 10 000 p–p collision initial state events at $\sqrt{s}=13\,\text{TeV}$. I map all the initiated partons inside the range $\eta<1$ into the transverse plane at each event. After that, I smear the partons with a 2D Gaussian distribution with width $\sigma_{g}$ weighted by the energy of each parton. After calculating $r_{\text{rms}}$ and $\epsilon_{2}$ at each event, I compute $\sigma_{r}^{(\text{AMPT})}$ and $\Gamma_{2}^{\epsilon(\text{AMPT})}$. The results are tabulated in Table 1 for different choices of $\sigma_{g}$.

In the fourth column, I have presented the kurtosis of elliptic flow fluctuation,

[TABLE]

from AMPT+GubsHyd at $\langle N_{\text{ch}}\rangle\approx 150$. I compute this quantity (which is Eq. (17) for $v_{2}$ fluctuation) because the term $\chi\sigma_{\epsilon}$ is canceled out from the numerator and denominator, and at the same time its sign is the same as $c_{2}\{4\}$. A similar quantity called the effective number of sources, $N_{s}$, has been measured for p–p collisions by ATLAS Aaboud:2017blb and can be translated to kurtosis by $\Gamma_{2}^{v}=-4/(N_{s}+3)$. At charge multiplicity $\langle N_{\text{ch}}\rangle\approx 150$, one finds the following range for the ATLAS measurement:

[TABLE]

Referring to Table 1, one observes that the negative sign for kurtosis is more compatible with $\sigma_{g}=0.1\,\text{[fm]}$. A p–p collision hydrodynamic study uses the HIJING+iEBE-VISHNU model in which a positive sign for $c_{2}\{4\}$ has been obtained Zhao:2017rgg . In addition to examining two values $\sigma_{g}=0.2,0.4\,[\text{fm}]$ for the smearing width, a range of other initial state and hydrodynamic parameters is tested in this paper. For all cases, the simulation outcome for $c_{2}\{4\}$ (and consequently $\Gamma_{2}^{v}$) has a positive sign. The GubsHyd model with the parameters mentioned in Eq. (20) returns the value $\Gamma_{2}^{v(\text{Gubs})}=-0.41$ which is the most accurate result compared to AMPT+GubsHyd and HIJING+iEBE-VISHNU.

By comparing the content of Table 1 with Eq. (20), one observes tension between AMPT initial state fluctuation properties with those extracted from the measurement directly by employing GubsHyd (Eqs. (19)). The rms radius fluctuation is more compatible with smearing width $\sigma_{g}=0.4\,[\text{fm}]$ while kurtosis is more compatible with $\sigma_{g}=0.1\,[\text{fm}]$.

The initial state model T${}_{\text{R}}$ENTo with nucleon substructure successfully describes the initial state of large and small collision systems. I use this model with parameters calibrated by Bayesian analysis Moreland:2018gsh . For the calibration, experimental measurements of Pb–Pb and p–Pb collisions at $\sqrt{s_{\text{NN}}}=5.02\;$TeV have been used. A full description of the parameters can be found in Ref. Moreland:2018gsh . Here, I focus on two parameters that impact the initial state fluctuation: $n_{c}$, the number of constituents inside the nucleon, and $\chi_{\text{struct}}$ which controls the width of the constituents ($\chi_{\text{struct}}=0$ corresponds to a nucleon with well-separated hot spots and $\chi_{\text{struct}}=1$ corresponds to a nucleon as a single Gaussian blob). The calibrated values (maximum a posteriori) of these parameters are $n_{c}=6.3$ and $\chi_{\text{struct}}=0.35$. I generate 100 000 p–p collision initial state events and compute kurtosis of ellipticity fluctuation and rms radius fluctuation width:

[TABLE]

The result shows tension between T${}_{\text{R}}$ENTo initial state fluctuation and Eq. (20). I also compute the hydrodynamic response to the T${}_{\text{R}}$ENTo initial state by using GubsHyd. At $\langle N_{\text{ch}}\rangle\approx 150$, T${}_{\text{R}}$ENTo+GubsHyd leads to \Gamma_{2}^{v(\tt\text{T{}{\text{R}}EN}+Gubs)}=0.26, meaning $c_{2}\{4\}$ is positive. I checked all combinations of $n_{c}=2,4,6,8$ and $\chi_{\text{struct}}=0.1,0.3,0.5,0.7$. The largest negative kurtosis for $\epsilon_{2}$ fluctuation is obtained by $n_{c}=8$ and $\chi_{\text{struct}}=0.1$ as \Gamma_{2}^{\epsilon(\tt\text{T{}{\text{R}}EN})}\approx-0.18. For rms radius fluctuation width, one finds \sigma_{r}^{\text{T{}{\text{R}}EN}}\approx 0.37\,[fm] in this choice of parameters. Even in this case, the kurtosis of elliptic flow fluctuation is positive, \Gamma_{2}^{v(\tt\text{T{}{\text{R}}EN}+Gubs)}\approx 0.18.

A more realistic hydrodynamic simulation that includes more stages (namely, pre-equilibrium and afterburner) is needed for a rigorous conclusion. Given that the estimations based on AMPT and T${}_{\text{R}}$ENTo could not predict a correct sign for $c_{2}\{4\}$ irrespective of its actual value, one might deduce that there is a missing piece in modeling the initial states or in modeling the collective evolution in p–p collisions. I will elaborate on this conclusion in the next section.

V Conclusion

In this paper, I introduced GubsHyd, a simple hydrodynamic “model” based on the Gubser solution. I showed that there is a lower bound for the size of a hydrodynamized system. I found that the flow measurements of p–p collisions are inside the domain of validity. I checked the validity of the presented simple model by comparing MC-Glauber+GubsHyd and MC-Glauber+VISH2+1 predictions in Pb–Pb collisions. I discussed the limitations of GubsHyd that should be taken into account in the conclusions. I found that three normalization/moments of the initial energy density, namely total energy density in the transverse direction, rms radius, and ellipticity (see Eqs. (2)), approximate the initial energy density and lead to a reasonable agreement between GubsHyd and iEBE-VISHNU for elliptic flow predictions in a wide range of centrality classes. The model could explain $v_{2}\{2\}$ and $v_{2}\{4\}$ measured by CMS and ATLAS collaborations for high multiplicity p–p collisions. Moreover, for systems close to the smallest QCD droplet, the model predicts an experimental signal in $v_{2}\{2\}$ and $v_{2}\{4\}$ measurements at lower multiplicities where in-plane/out-of-plane hydrodynamic response transition happens.

Employing the simplicity of the Gubser solution, I found a smooth function for a hydrodynamic response to the initial ellipticity. This response function is given in terms of final state multiplicity and rms radius of the initial state $k_{2}(n_{\text{tot}},r_{\text{rms}})$ (see Fig. 3). Given that $\epsilon_{2}$ and $r_{\text{rms}}$ play the most significant roles in the final value of elliptic flow, I modeled ellipticity and rms radius fluctuations. Then, using the smooth function $k_{2}(n_{\text{tot}},r_{\text{rms}})$, I explicitly obtained elliptic flow fluctuation in terms of initial state fluctuation properties (Eqs. (19)). Since I do not use any specific model for the initial state energy density, I could extract the fluctuation properties from the experimental measurements. By comparing GubsHyd predictions for two- and four-particle correlations with those measured for p–p collisions, I found that the kurtosis of ellipticity fluctuation should be -3/4 approximately (see Eq. (20)). I examined two initial state models, AMPT and T${}_{\text{R}}$ENTo, in a wide range of parameters. None of these models could reproduce such a significant negative value for the kurtosis and, at the same time, a reasonable value for rms radius fluctuation width. This observation could explain the unsuccessful hydrodynamic predictions for $c_{2}\{4\}$ in p–p collisions in recent studies.

Regarding the sign of $c_{2}\{4\}$, a more realistic hydrodynamic simulation is essential for the future. However, considering the validation of GubsHyd with iEBE-VISHNU in Sec. III and appendix B, I would expect that the conclusions based on GubsHyd should be close to a conventional causal hydrodynamic study. One notes that the effect of hydrodynamization surface becomes important at lower multiplicities. The four-particle correlation, $c_{2}\{4\}$, is measured at high enough multiplicities that one can ignore the impact of the hydrodynamization surface. These remarks bring one to conclude the following scenarios: First, the anisotropic flow in p–p collisions can be explained by conventional hydrodynamics, but there is a missing piece in the initial state models that prevents the model from producing a large negative kurtosis for ellipticity fluctuation. Second, conventional hydrodynamics alone cannot explain p–p collision observations while the initial states have the correct fluctuating properties. One candidate for the second scenario could be the presence of hydrodynamic fluctuations Kapusta:2011gt . In particular, in the context of Gubser flow, it has been demonstrated that hydrodynamic fluctuations have non-negligible effects on flow measurements of p–p collisions Yan:2015lfa . This scenario is worth a more careful study in the future.

ACKNOWLEDGMENT

I thank A. Bilandzic for comments and discussions. I also thank C. Mordasini and M. Lesch for comments. I thank U. A. Wiedemann, A. Kurkela, and W. van der Schee for discussions. This project has received funding from the European Research Council (ERC) under the European UnionsHorizon 2020 research and innovation program (Grant No. 759257).

Appendix A A simple hydrodynamic model based on the Gubser solution (GubsHyd)

The initiation and hydrodynamic evolution in the GubsHyd model are based on the Gubser solution. Although the solution is briefly mentioned in the main text, I present it with slightly more details as it is necessary for the rest of the appendix.

Instead of $\mathbb{R}^{1,3}$ space, one is able to employ the conformal symmetry to represent a generic field $\Phi$ in $d\hat{s}^{2}=\Omega^{-2}ds^{2}$ space with $\Omega=\tau$. I show $\Phi$ in the rescaled space $d\hat{s}^{2}$ with a hat. Depending on the mass dimension and the type of the tensor $\Phi$, the appropriate scaling should be employed as well. For instance, a scalar quantity $\phi$ with mass dimension $\Delta$ is written in $d\hat{s}^{2}$ space as $\hat{\phi}=\tau^{\Delta}\phi$. In this space, the Gubser solution with an elliptic perturbation on top of that is given by Gubser:2010ze ; Gubser:2010ui ; Staig:2011wj

[TABLE]

where

[TABLE]

and $Y_{l,m}\equiv Y_{l,m}(\theta,\phi)$ is the spherical harmonic function. The hydrodynamic equations are given by $\partial_{\mu}T^{\mu\nu}=0$ where

[TABLE]

is the ideal stress-energy tensor and $\eta^{\mu\nu}=\text{diag}(-1,1,1,1)$ is the metric. Substituting the ansatz in Eq. (25) into hydrodynamic equations and using conformal e.o.s. $\epsilon=3p$, one obtains the following system of differential equations for $\delta_{2}(\rho)$ and $\tilde{\delta}_{2}(\rho)$:

[TABLE]

Initial isotropic in fluid velocity assumption fixes the initial condition as

[TABLE]

This solution can be translated to hyperbolic-cylindrical coordinates via

[TABLE]

The connection between $(q,\hat{\epsilon}_{0},\lambda,\rho_{\text{hyd}})$ and standard parameters $(r_{\text{rms}},\epsilon_{\text{tot}},\epsilon_{2},\tau_{\text{hyd}})$ has been discussed in the main text. For the next step, I compute the hadronization of the above solution.

I focus on events with size $r_{\text{crit}}<r_{\text{rms}}<R_{\text{crit}}$ in which a hydrodynamized matter, “QGP”, is formed at least in a small region in the energy density (core) while at the tail of the energy density (corona), the system is in the hadronic phase. The spectrum of such a system can be written as ($d\boldsymbol{p}\equiv d^{3}p/E$)

[TABLE]

I assume the core part evolves with hydrodynamic equations, while free streaming starts immediately after the initiation for the corona part.

I employ the Cooper-Frye prescription [see Eq. (5)] to obtain the particle distribution ($g=2(N_{c}^{2}-1)=16$ has been chosen). In this paper, I have considered Boltzmann equilibrium distribution instead of a more sophisticated Bose-Einstein/Fermi-Dirac near-equilibrium distribution.

For the core part, the particles are emitted after the hydrodynamic evolution from freeze-out surface $\Sigma^{\mu}=(\rho,\theta_{\text{fo}}(\rho,\phi),\phi,\eta)$ with surface element

[TABLE]

The functionality of $\cos\theta_{\text{fo}}(\rho,\phi)$ for the Gubser background solution ($\lambda=0$) has been mentioned in Eq. (6). To obtain anisotropy in the final particle distribution, however, one needs to keep the elliptic perturbation. This modifies the Eq. (6) to the following:

[TABLE]

where

[TABLE]

and $\theta_{b}(\rho)=\cos\theta_{\text{fo}}(\rho,\phi)|_{\lambda=0}$.

I start with computing the final multiplicity. Because the Gubser solution has a more straightforward analytical form in rescaled space $d\hat{s}^{2}$, I first calculate Eq. (5) in the rescaled space and then I transform $dN_{\text{core}}/d\boldsymbol{\hat{p}}$ into $dN_{\text{core}}/d\boldsymbol{p}$, using the fact that the latter has mass dimension -2 (recall $d\boldsymbol{p}\equiv d^{3}p/E$):

[TABLE]

One notes that the four-momentum $p^{\mu}$ ($p_{\mu}p^{\mu}=-m^{2}$) in the hyperbolic-cylindrical coordinate reads as

[TABLE]

where $(p_{T},\phi_{p})$ is the transverse momentum in the polar coordinate, $m_{T}$ is the transverse mass, and $y_{p}$ is rapidity. Subsequently, one obtains the momentum in $d\hat{s}^{2}$ via

[TABLE]

where $\hat{x}^{\mu}=(\rho,\theta,\phi,\eta)$ and $x^{\mu}=(\tau,r,\phi,\eta)$. Keeping terms up to linear order in $\lambda$, the particle distribution for the core part is written as

[TABLE]

where

[TABLE]

The extra $\tau_{\text{fo}}=\tau(\rho,\theta_{\text{fo}})$ in the measure is coming from Eq. (35). The relation $\hat{T}_{\text{fo}}=\tau_{\text{fo}}T_{\text{fo}}$ has also been used.

To calculate the multiplicity, I concentrate on $\lambda=0$ (background solution). In this case, the integration on $\eta$ and $\phi$ in Eq. (38) can be performed trivially. The result reads as

[TABLE]

where

[TABLE]

The particle distribution of the corona part can also be obtained by the Cooper-Frye formula, where the freeze-out surface is coincident with the hydrodynamization surface, $\Sigma^{\mu}=(\rho_{\text{hyd}},\theta,\phi,\eta)$:

[TABLE]

where $\hat{T}(\rho_{\text{hyd}})=(\hat{\epsilon}_{b}/C_{0})^{1/4}/\cosh^{2/3}\rho_{\text{hyd}}$ is obtained from (25) at $\lambda=0$ and e.o.s. $\epsilon=C_{0}T^{4}$. Similar to the core part, one can trivially perform the integration on $\phi$ and $\eta$:

[TABLE]

where

[TABLE]

The integrals in Eqs. (40) and (43) can be computed numerically. Referring to Eq. (31), one obtains the total multiplicity, $dN_{\text{QGP}}/d\boldsymbol{p}|_{\lambda=0}$, by adding the outcome of two integrals. I have prepared the integrals in a Mathematica package where the $p_{T}$ dependent and $p_{T}$ integrated multiplicity in the unit rapidity are available:

[TABLE]

The total multiplicity in the unit rapidity $n_{\text{tot}}(r_{\text{rms}},\epsilon_{\text{tot}})$ is obtained from the second function with $p_{T}$ range from zero to infinity. One can also change the parameters of final particle distribution from $y_{p}$ to space-time pseudorapidity via

[TABLE]

Now by using the above and by inverting $n_{\text{tot}}(r_{\text{rms}},\epsilon_{\text{tot}})$ to find $\epsilon_{\text{tot}}(n_{\text{tot}},r_{\text{rms}})$, one obtains $dN_{\text{QGP}}/p_{T}dp_{T}d\eta|_{\lambda=0}$ and $N_{\text{tot}}|_{\lambda=0}$ as a function of total multiplicity in the unit rapidity, rms radius, transverse momentum, and space-time pseudorapidity:

[TABLE]

Having found the energy dependence of multiplicity, one can translate the lower and upper bound of the system size into multiplicity. The total energy density for a system at critical size can be obtained from Eq. (7),

[TABLE]

Substituting $\epsilon_{\text{crit}}$ into $n_{\text{tot}}(r_{\text{rms}},\epsilon_{\text{tot}})$, one obtains a critical value for the multiplicity. Another critical value can be found by inverting $R_{\text{crit}}(\epsilon_{\text{tot}})$ (see Eq. (8)) and substituting the result into $n_{\text{tot}}(r_{\text{rms}},\epsilon_{\text{tot}})$. Both of these critical multiplicities are available in the package,

[TABLE]

The above functions are shown with thick black curves in Fig. (3) for $0.3<p_{T}<3\,\text{GeV}$. In massless limit $m=0$, both $n_{\text{tot}}(r_{\text{rms}},\epsilon_{\text{tot}})$ and $n_{\text{crit}}$ (I assume a full range $p_{T}$ integration) can be obtained analytically,

[TABLE]

As can be seen in the massless limit, $n_{\text{crit}}$ has no $r_{\text{rms}}$ dependence.

Now, one can compute $k_{2}(r_{\text{rms}},n_{\text{tot}})=\left.(\partial v_{2}/\partial\epsilon_{2})\right|_{\epsilon_{2}=0}$. The final particle distribution is obtained via

[TABLE]

The first term in the right-hand side is already computed. In the second term, only the core part contributes to momentum anisotropy at the final state. This is because I have assumed free streaming for the corona part. Therefore, the elliptic flow can be obtained by

[TABLE]

One can write the linearized $dN_{\text{core}}/d\boldsymbol{p}$ in Eq. (38) with respect to $\lambda$ as follows:

[TABLE]

where

[TABLE]

and

[TABLE]

For a given function $f(x,t)$, I employ the following calculus identity,

[TABLE]

By using the above identity, the derivative of Eq. (53) at $\lambda=0$ is given by

[TABLE]

where the prime in $\rho^{\prime}_{\text{hyd}}(\lambda)$ is the derivative with respect to $\lambda$. One notes that, by definition, $\cos\theta_{\text{fo}}(\rho_{\text{max}})=1$, which immediately leads to $\sin\theta_{\text{fo}}(\rho_{\text{max}})=0$. Given that there is a term containing $\sin\theta_{\text{fo}}$ in $\mathcal{Q}_{0}$, the boundary term $\mathcal{Q}_{0}e^{\chi_{0}}|_{\rho_{\text{max}}(0)}$ in Eq. (57) is identically vanishing. The second boundary term, however, is not zero. To find $\rho^{\prime}_{\text{hyd}}(0)$, one needs to compute $\rho_{\text{hyd}}$ at finite $\lambda$. The quantity $\hat{\epsilon}_{0}$ up to the linear term in $\lambda$ is given by

[TABLE]

Noting $\epsilon=\hat{\epsilon}/\tau^{4}$ and substituting Eq. (58) into Eq. (25a), one obtains

[TABLE]

where in the above $\theta_{\text{fo}}\equiv\theta_{\text{fo}}(\rho_{\text{hyd}}(0),\phi)$. One can compute $\rho^{\prime}_{\text{hyd}}(0)$ from Eq. (59).

In Eq. (57), the integration over $\phi$ and $\eta$ is tedious but analytically doable, while the integration over $\rho$ needs to be done numerically. Symbolically, one can write Eq. (57) as follows:

[TABLE]

where $\mathcal{V}_{0}(p_{T})$ and $\mathcal{V}_{2}(p_{T})$ are two complicated functions depending on the standard quantities including a numerical integration over $\rho$. Using Eq. (52), one finds

[TABLE]

where I have substituted $\lambda=\left(\sqrt{5\pi}/3\right)\epsilon_{2}$. This function and its $p_{T}$ integrated version are available in the Mathematica package,

[TABLE]

Also, by substituting $\epsilon_{\text{tot}}(n_{\text{tot}},r_{\text{rms}})$ into $k_{2}(\epsilon_{\text{tot}},r_{\text{rms}})$, I obtain hydrodynamic response in terms of rms radius and multiplicity in the unit rapidity:

[TABLE]

The second function in the above within range $0.3<p_{T}<3\,\text{GeV/c}$ has been employed to plot Fig. (3).

Appendix B Validation of GubsHyd via iEBE-VISHNU

The dissipative effects in both hydrodynamic evolution and freeze-out have been ignored in the GubsHyd model. In the main text, I assume that the contribution of dissipative effects on integrated $v_{2}\{2\}$ and $v_{2}\{4\}$ can be corrected via an overall factor $\chi$ (see Eq. (12)). I have examined this assumption by comparing the prediction of MC-Glauber+GubsHyd and MC-Glauber+VISH2+1 (iEBE-VISHNU) with nonvanishing shear viscosity over entropy density for integrated $v_{2}\{2\}$ and $v_{2}\{4\}$ of charged pions (see Fig. 4). In this appendix, I extend the comparison to the following observables: $p_{T}$ spectrum, differential $v_{2}\{2\}$, and $v_{2}\{4\}$ and centrality dependence of the multiplicity for charged pions. These comparisons lead to a clearer picture of the accuracy of GubsHyd and the limitations that should be considered in the conclusions (see Ref. Dusling:2007gi for a comprehensive study about the role of the dissipative effects on anisotropic flow and $p_{T}$ spectrum).

The contribution of dissipative effects on elliptic flow can be explained as follows: In the course of evolution, compared to an ideal fluid, the initial elliptic shape of the energy density becomes more round at the presence of the dissipative effects. Consequently, $v_{2}$ fluctuates in a narrower range event by event. This argument explains that $\chi$ in Eq. (12) should be smaller than unity to get a reasonable agreement between GubsHyd and iEBE-VISHNU in Fig. 4. One notes that the dissipative effects are dynamical and cannot be fully corrected via a simple scaling factor. In Fig. 6, the $p_{T}$ dependence of $v_{2}\{2\}$ and $v_{2}\{4\}$ is shown where parameters are tuned as those mentioned in Sec. III. As seen, the GubsHyd prediction starts deviating from iEBE-VISHNU prediction above $p_{T}\sim 1$GeV/c, meaning the dissipative effects have more influence on the harder particles. Since the particles with $p_{T}\lesssim 1$GeV/c are more abandoned, the part with $p_{T}\gtrsim 1$GeV/c is less important in integrated $v_{2}\{2k\}$ observables. As a result, a constant $\chi$ leads to a reasonable approximation for integrated $v_{2}\{2\}$ and $v_{2}\{4\}$.

Now, I discuss the $p_{T}$ spectrum and multiplicity of charged pions. The dissipative effects lead to entropy production and consequently increase the number of particles in the final state. As a result, one expects the following relation:

[TABLE]

where $\chi^{\prime}$ should be greater than unity. The numerical value of $\chi^{\prime}$ has no impact on the observables I discuss in the main text, but I analyze it here for completeness. Although the dissipation effects act dynamically on multiplicity and cannot be fully corrected by a constant $\chi^{\prime}$, it is still a good approximation to assume it is a constant in a wide range of centralities, similar to integrated $v_{2}\{2\}$ and $v_{2}\{4\}$. By demanding that the average of the number of charged pions in the unit rapidity, $\langle dN_{\pi^{\pm}}/dy_{p}\rangle$, is the same in the centrality range 0$-$80% for GubsHyd and iEBE-VISHNU, one finds $\chi^{\prime}\approx 1.4$. Using this number, I have plotted the $p_{T}$ spectrum of charged pions in Fig. 7 (left). As seen from the figure, the $p_{T}$ spectrum of iEBE-VISHNU has a smaller slope compared to GubsHyd. The reason is that, similar to differential $v_{2}$, the dissipative effects are more significant at higher $p_{T}$. As a result, there is more entropy production and more particles at higher $p_{T}$. The centrality dependence of charged pion multiplicity is depicted in Fig. 7 (right). One observes that the overall trends of GubsHyd and iEBE-VISHNU are in a reasonable agreement even though the prediction from iEBE-VISHNU shows a steeper trend.