Searching for small droplets of hydrodynamic fluid in proton--proton collisions at the LHC

In this paper, we investigate the hydrodynamic collectivity in high-multiplicity events of proton-proton collisions at $\sqrt{s}=$ 13 TeV, using iEBE-VISHNU hybrid model with three different initial conditions, namely, HIJING, super-MC and TRENTo. With properly tuned parameters, hydrodynamic simulations with each initial model give reasonable descriptions of the measured two-particle correlations, including the integrated and $p_{\rm T}$-differential flow for all charged and identified hadrons. However, the hydrodynamic simulations fail to describe the negative value of the four-particle cumulant $c_2^v\{4\}$ as measured in experiments. Further investigations show that the non-linear response between the elliptic flow $v_2$ and the initial eccentricity $\varepsilon_2$ becomes significant in the small p-p systems. This leads to a large deviation from linear eccentricity scaling and generates additional flow fluctuations, which results in a positive $c_2^v\{4\}$ even with a negative $c_2^\varepsilon\{4\}$ from the initial state. We also presented the first hydrodynamic calculations of multi-particle mixed harmonic azimuthal correlations in p-p collisions, such as normalized asymmetric cumulant $nac_n\{3\}$, normalized Symmetric-Cumulant, $nsc_{2,3}\{4\}$ and $nsc_{2,4}\{4\}$. Although many qualitative features are reproduced by the hydrodynamic simulations with chosen parameters, the measured negative $nsc_{2,3}\{4\}$ cannot be reproduced. The failure of the description of negative $c_2\{4\}$ and $nsc_{2,3}\{4\}$ triggers the question on whether hydrodynamics with a fundamentally new initial state model could solve this puzzle, or hydrodynamics itself might not be the appreciated mechanism of the observed collectivity in p-p collisions at the LHC.

These observed flow-like signals in the small systems can be quantitatively or semi-quantitatively described by * You Zhou: you.zhou@cern.ch † Huichao Song: huichaosong@pku.edu.cn hydrodynamic calculations [38][39][40][41][42][43][44][45][46][47][48][49][50][51], which translate initial spatial anisotropies into final momentum anisotropies of produced hadrons with the collective expansion of the bulk matter. Besides, other model calculations based on final state interactions, such as transport models [52][53][54][55][56], hadronic rescatterings [57,58], a string rope and shoving mechanism [59] have also been performed to study the collective behavior of small systems. Alternatively, the color glass condensate (CGC) or IP-Plasma focused on initial state effects [60][61][62][63][64][65][66][67][68][69] can also qualitatively reproduce many features of collectivity. The origin of the observed collective behavior in the small systems is still under intense debate. Recently, the model calculations in [70] showed that the quark coalescence procedure is necessary to reproduce the number of constituent quark scaling of v 2 at intermediate p T in high-multiplicity p-Pb collisions at √ s =13 TeV, which demonstrate the importance of the partonic degrees of freedom and possible formation of QGP in the small p-Pb systems.
[ 30], have been applied to remove the non-flow contaminations for the extracted flow harmonics. Meanwhile, multi-particle cumulants have been systematically measured, which provide more insights for the collective phenomenon in high-multiplicity p-p collisions. Compared to the two-particle correlations, multi-particle cumulants, by construction, have the advantage of suppressing shortrange two-particle correlations [71,72]. Besides, twoand three-subevent methods have been implemented to further suppress the remaining non-flow contaminations, which are also much less sensitive to the multiplicity fluctuations compared to the standard method [72][73][74]. It was found that c 2 {4} turns to negative value in highmultiplicity events of p-p collisions, which gives the real value of the flow coefficients v 2 {4} through the relation c 2 {4} = −v 2 {4} 4 and strongly indicates the existence of anisotropic flow in the small p-p systems [73,74]. Furthermore, ALICE [32], ATLAS [74] and CMS [75] have measured the correlations between different flow harmonics v n and v m via three-or four-particle cumulants in p-p collisions, which shows negative correlations between v 2 and v 3 and positive correlations between v 2 and v 4 with similar relative correlation strengths as measured in p-Pb and Pb-Pb systems. It is thus on-time and important to investigate these collective flow signatures by hydrodynamic models and to discuss whether hydrodynamic calculations could describe two-and multi-particle cumulants simultaneously in the small p-p systems created at the LHC.
In our previous work [76], we found that, with properly tuned parameters, hydrodynamic simulations with HIJING initial conditions can nicely describe the twoparticle correlations in p-p collisions at √ s =13 TeV, including the integrated v 2 {2}, differential elliptic flow v 2 (p T ) for all charged hadrons and for identified particles (K 0 S and Λ). However, the measured negative c 2 {4}, which has been usually interpreted as evidence of hydrodynamic flow, could not be reproduced by our hydrodynamic calculations which showed a positive value of c 2 {4}. It is still unknown if the wrong sign of c 2 {4} is due to the incorrect initial conditions from HIJING or due to the application of the hydrodynamic model to p-p collisions.
To address these questions, in this paper, we implement three different initial conditions, called HIJING [51], super-MC [77] and TRENTo [78], to the iEBE-VISHNU hybrid model simulations to study various flow observables in p-p collisions, especially on the four-particle cumulant c 2 {4} and mixed harmonic three-and fourparticle azimuthal correlations. To better understand the non-linear hydrodynamic evolution in the small systems, we also investigate the response between the initial ε 2 and final v 2 . In addition, we study the effects of pre-equilibrium dynamics in the p-p collision by including the free-streaming evolution before the hydrodynamic simulations.
This paper is organized as follows: in the next section, we will give an introduction of the iEBE-VISHNU hy-drodynamic model and the initial conditions of HIJING, super-MC and TRENTo and explain the setups. Section III presents the model calculations, the comparison to the experimental data, and the related discussion. Section IV gives a brief summary of this paper.

II. THE MODEL AND SET-UPS
A. iEBE-VISHNU hybrid model iEBE-VISHNU [79] is an event-by-event version of hybrid model VISHNU [80] that combines 2+1D viscous hydrodynamics VISH2+1 [81,82] to describe the QGP expansion with a hadron cascades model UrQMD [83,84] to simulate the evolution of hadronic matter. Based on the Israel-Stewart formalism, VISH2+1 solves the transport equations for the energy-momentum tensor T µν and shear stress tensor π µν with a state-of-the-art equation of state (EoS) s95-PCE [85,86] as an input to simulate viscous fluid expansion of the hot QCD matter with longitudinal boost-invariance. For simplicity, we neglect the bulk viscosity, net baryon density, and heat conductivity and assume a constant specific shear viscosity η/s. The hydrodynamic evolution matches the hadron cascade simulations at a switching temperature T sw , where various hadrons are emitted from the switching hypersurface for the succeeding UrQMD evolution.
To systematically investigate the hydrodynamic collectivity of p-p system and its dependence on the initial condition models, we implement three different initial condition model, namely, modified HIJING [51], super-MC [77] and TRENTo [78]. In general, these three initial conditions neglect the pre-equilibrium dynamics and set the initial flow velocity and the shear-stress tensor to be zero for the succeeding hydrodynamic simulations, which are also the default settings of our calculations. In this paper, for one parameter set of TRENTo initial condition, we prepared a version including the free-streaming evolution before hydrodynamics to study the pre-equilibrium effects. Below is a brief description of these three initial condition models.

B. HIJING initial condition
In HIJING [87][88][89], the radial density of the colliding protons is the Woods-Saxon shapes, and the produced jet pairs and excited nucleus are treated as independent strings, where the hard jet productions are calculated by pQCD, and the soft interactions are treated as gluon exchange within Lund string model. For the HIJING initial condition developed in Ref. [51], it assumes that the mother strings that break into independent partons quickly form several hot spots for the succeeding hydrodynamic evolution. The center positions of these mother strings (x c , y c ) are sampled by the Woods-Saxon distribution, and the positions of the produced partons (x i , y i ) within the strings are sampled with a Gaussian distribution with a width σ R : exp[− (xi−xc) 2 +(yi−yc) 2 The initial energy density profiles in the transverse plane for the 2+1D hydrodynamic evolution are constructed from the energy depositions of emitted partons, together with an additional Gaussian smearing [51] e(x, y) = K where σ 0 is the Gaussian smearing factor, p i is the momentum of the produced parton i, and K is an additional normalization factor. Here, we neglect the initial flow U 0 and only consider the partons within the mid-rapidity |η| < 1 (for related details, please also refer to Ref. [51]).

C. super-MC initial condition
For p-p collisions, super-MC model with sub-nucleonic fluctuations [77] assumes the colliding protons consist of three valence quarks, and the collisions between valence quarks depose a fraction of the kinetic energy of the colliding systems into the initial energy of the newly formed matter, which fluctuates from event to event. Following Ref. [77], the initial entropy density of the produced matter is modeled as where γ (i) k (i = 1, 2, 3) are the random weighting factors that used to fit the multiplicity distributions in p-p collisions, r i (i = 1, 2, 3) is the position of three valence quarks which is distributed according to a Gaussian probability distribution. σ g is a factor to describe the shape of quark density distribution together with a consideration of low-x gluon contributions [77].

D. TRENTo initial condition
TRENTo is a parameterized initial condition model, which generates the initial entropy density via the reduced thickness function [90,91]: whereT (x, y) is the modified participant thickness function, s 0 is a normalization factor, and p is a tunable parameter which makes TRENTo model effectively interpolates among different entropy deposition schemes, such as KLN, EKRT, WN, etc. [78,90,91]. For proton-proton collisions, TRENTo is modified with the sub-nucleonic structure [78] so thatT (x, y) is written where n c is the number of independent constituents in a proton, γ i (i = 1, 2, ..., n c ) is a random weighting factor with the unit mean and variance 1/k, x i (i = 1, 2, ..., n c ) are the positions of constituents, b is the impact parameter, and ρ c is the density of constituents written in a Gaussian form: 2v 2 ), and v is a tunable effective width of nucleons.
TRENTo initial condition with free streaming: For TRENTo initial condition, we also construct another type of the initial condition with free-streaming to include the effects of pre-equilibrium dynamics before hydrodynamic evolution. Following Refs. [78,92,93], we assume the particle density of non-interacting massless particles at the very beginning is proportional to entropy density described by Eq. (3), and then free streaming these massless partons till the proper time τ 0 to obtain the boost-invariant energy-momentum tensor T µν (x, y, τ 0 ). After that, we implement the following Landau matching condition to obtain the initial energy density e(x, y, τ 0 ) and fluid velocity u µ (x, y, τ 0 ): and the initial shear stress tensor and bulk pressure can be calculated with:  and TRENTo (c) initial conditions. The CMS and ATLAS data are taken from Refs. [30,75] and Refs. [31,74], respectively.
with the spatial projector being ∆ µν = g µν − u µ u ν , together with an equation of state of P = 1 3 e for the massless ideal gas at the initial state.
For p-p collisions at √ s =13 TeV, we implemented several sets of parameters for each of these three or four different initial conditions. These parameters are roughly tuned to approximately fit the p T -spectra [94] and v 2 {2} [30,31,74,75] measured in experiments. Note that these data are not enough to fully constrain the free parameters in hydrodynamic simulations. Since this paper is aimed to investigate the sign of c 2 {4}, mixed harmonic three-and four-particle azimuthal correlations, and the effects of non-linear evolution rather than make quantitative descriptions and prediction for p-p collisions, we chose three or four sets of parameters for each initial condition, as listed in Tables I, II and III. For TRENTo initial condition with the parameter set Para-I, we consider two cases with and without free-streaming as described above.

A. 2-particle cumulant
With various sets of parameters for these three initial conditions, HIJING, super-MC and TRENTo, as listed in Tables I, II and III, we calculate the integrated v n {2} (n = 2, 3 and 4) as a function of multiplicity for pp collisions at √ s = 13 TeV, using iEBE-VISHNU together with an application of the two-subevent method with the pseudorapidity gap |∆η| > 0, kinematic cuts 0.3 < p T < 3.0 GeV/c and |η| < 2.4. To eliminate the effects of multiplicity fluctuations, we implement the same method as used in experimental analysis and in our early paper [31,51], which first obtain the 2-and 4-particle cumulants within the multiplicity class with the number of charged hadrons N sel ch with 0.3 < p T < 3.0 GeV/c and |η| < 2.4, and then map it to the number of charged hadrons N ch with 0.4 < p T GeV/c and |η| < 2.4 to compare with the experimental data. Fig. 1 presents the comparison between our hydrodynamic calculations and the experimental measurements from ATLAS [30,75] and CMS [31,74]. It shows that hydrodynamic simulations with these three different initial conditions can generally reproduce the multiplicity dependence of the integrated v 2 {2} as we could expect from tuning the related parameters. Note that these four sets of parameters, Para-I-IV, in HIJING initial condition, are the same as we used in Ref. [51], which are tuned to fit v 2 {2} data obtained from the "peripheral subtraction" method (Para-I-III) and from the "template fit" method (Para-IV), respectively. For super-MC and TRENTo initial conditions, we choose one set of parameters (Para-III for super-MC and Para-II for TRENTo) to describe the v 2 {2} data with "peripheral subtraction" method, and the other parameter sets to approximately describe the data with "template fit" method. In general, hydrodynamic calculations approximately describe v 4 {2} from CMS and ATLAS, but tend to overestimate the measured v 3 {2} with both "peripheral subtraction" and "template fit" methods, especially for the ones obtained with TRENTo initial conditions. On the other hand, v 3 {2} data from "peripheral subtraction" and "template fit" methods also largely deviate from each other, and it is still under debate on which method gives a better non-flow subtraction for the odd flow harmonics [74].
For the parameter set of Para-I of TRENTo initial condition, we also include the pre-equilibrium evolution with an infinitely weak coupling limit (dashed red line), √ s = 13 TeV, calculated by iEBE-VISHNU with HIJING, super-MC and TRENTo initial conditions. The CMS and ATLAS data are taken from Refs. [30] and [29], respectively. which free-streams the initial state to proper time τ 0 before instantaneously switching to hydrodynamic simulations [78,93]. It shows that such pre-equilibrium dynamics not only affects the magnitude of v 2 {2} but also affects its dependence on the multiplicity, which seems excluded by the experimental data.
From the hydrodynamic calculations shown in Fig. 1, it is clear that the flow coefficients of v 2 , v 3 and v 4 in p-p collisions could provide certain constraints on the parameter settings for model calculations with various initial conditions. A simultaneous description of v 2 , v 3 and v 4 is one of essential steps to validate the applicability of hydrodynamic simulations in small systems.
In Fig. 2, we calculate differential elliptic flow v 2 (p T ) for all charged hadrons (a)-(c) and for K 0 S and Λ (d)-(f) for the multiplicity range 80 < N sel ch < 120 with the two-particle cumulant method with a pseudorapidity gap |∆η| > 0. iEBE-VISHNU with HIJING, super-MC or TRENTo initial conditions can roughly describe v 2 (p T ) for all charged hadrons measured from CMS and ATLAS with the "peripheral subtraction" or "tem-plate fit" method. More specifically, as mentioned previously in Ref. [51], the calculations of Para-I, II and III of HIJING initial condition, which are tuned for "peripheral subtraction", give a satisfactory description of the data. In contrast, Para-IV of HIJING initial condition tuned for "template fit" slightly overpredicts the data above 1.0 GeV/c. For super-MC and TRENTo initial conditions, hydrodynamic calculations can roughly describe v 2 (p T ) data for all charged hadrons.
The panels (d)-(f) of Fig. 2 present v 2 (p T ) for identified hadrons, which show clear v 2 mass ordering between K 0 S and Λ for both CMS measurements and our iEBE-VISHNU calculations. The hydrodynamic predictions with HIJING initial condition (Para-I and II) can nicely describe the data. However, the calculations with super-MC initial condition tend to overestimate v 2 (p T ) of K 0 s and Λ, and the calculations with TRENTo initial condition tend to underestimates the data of Λ. The mass splitting between K 0 S and Λ is more significant for calculations with TRENTo initial condition. Such larger mass splitting of v 2 indicates a stronger radial flow de- and TRENTo (c) initial conditions using standard cumulant method. The CMS data with standard cumulant method and the ATLAS data with three-subevent method are taken from Refs. [30] and [31], respectively. velopment during the hydrodynamic evolution. This is consistent with what we have seen (but not shown here) in the p T -spectra (the spectra obtained with TRENTo initial condition is harder than the others [95]), which can also provide certain constrains on the initial conditions.

B. 4-particle cumulant
In Fig. 3, we study the four-particle cumulants of the second harmonics, c 2 {4}, in high-multiplicity proton-proton collisions at √ s = 13 TeV. Although iEBE-VISHNU can roughly describe the measured v n {2} using these three initial conditions with the properly tuned parameters, the predicted c 2 {4} are always positive in the high-multiplicity region and fail to reproduce the negative c 2 {4} as measured in experiments. In Ref. [51] we have found that the positive c 2 {4} from hydrodynamic simulations with HIJING initial condition is not due to the effects of non-flow contributions or multiplicity fluctuations. We also demonstrated that the standard method, two-subevent method and three-subevent method almost give the same value of c 2 {4} for such flow-dominated systems. The panels (b) and (c) in Fig. 3 also show, for the two newly implemented super-MC and TRENTo initial conditions, iEBE-VISHNU still generates a positive c v 2 {4} even for these parameter sets associated with a negative c ε 2 {4} for 0-0.1% events in the initial states, as listed in Table IV. Note that recent MUSIC hydrodynamic simulations with IP-Glasma initial conditions also give positive c v 2 {4} for the entire multiplicity range in p-p collisions at √ s = 13 TeV [96]. We thus emphasize that hydro- With such findings, we then focus on the effects of non-linear hydrodynamic evolution on the four-particle cumulant c 2 {4}. Specifically, if the final v 2 has a linear response to the initial ε 2 , the scaled v 2 distributions P (v 2 / v 2 ) should overlap with the scaled ε 2 distributions P (ε 2 / ε 2 ), which is the case for central and midcentral Pb-Pb collisions [97]. If such a linear response holds in p-p collisions, the final state c v 2 {4} is expected to have the same sign as the initial state c ε 2 {4}. However, hydrodynamic simulations did not confirm such such expectation. As shown in Figs. 3 and 4, and Table IV, the negative initial c ε 2 {4} (e.g., Para-III of HIJING, Para-I of super-MC, and Para-I-III of TRENTo) still lead to a positive c v 2 {4} at final state after the hydrodynamic evolution. We also find that even though the c ε 2 {4} of Para-III of TRENTo initial conditions is more negative than that of Para-I in Table IV, the initial ε 2 distribution of Para-III of TRENTo initial condition is "wider" with larger mean value of ε 2 . The corresponding larger non-linear effects during the evolution lead to a larger positive value of c v 2 {4} than the one associated with Para-I. As demonstrated by Figs. 3 and 4 and also confirmed by additional calculations which are not shown in this paper, similar situations also happen for iEBE-VISHNU simulations with HIJING or super-MC initial conditions. To further understand the general "wrong sign" of c v 2 {4} from hydrodynamic simulations with various initial conditions, we study the correlation between initial eccentricity ε 2 and final elliptic flow v 2 . As shown in Fig. 5 (a), a clear deviation of elliptic flow from linear scaling is observed for ε 2 > 0.5 where the cubic term becomes significant, which is similar to the peripheral Pb-Pb collisions [98] 1 . Such non-negligible cubic response leads to the fact that the scaled distribution P (v 2 / v 2 ) and P (ε 2 / ε 2 ) does not overlap with each other as shown in Fig. 5 (b). It also introduces additional fluctuations of v 2 in the final states, which could even change the sign of c v 2 {4} and make the model calculations fail to reproduce the negative c v 2 {4} measured in experiments. It has been generally argued that two-and multiparticle cumulants have different sensitivities to the flow fluctuations, which is written as [71] v Here v n and σ v represent the flow and flow fluctuations. These equations are valid in the case of small flow fluctuations, which might not be applied in small systems like p-p collisions. However, considering the fact that hydrodynamic calculations could quantitatively describe the two-particle correlations but could not even produce the correct sign of four-particle cumulants, one can conclude that the current hydrodynamic calculations could not simultaneously describe both the anisotropic flow v n and the flow fluctuations σ v .
1 It requires significant amount of computational resources to obtain In Fig. 6, we further study the normalized threeand four-particle azimuthal correlations in highmultiplicity proton-proton collisions at √ s = 13 TeV. The three-particle asymmetric cumulant is defined as n v 2n cos 2n(Ψ n − Ψ 2n ) , which is sensitive to the correlations between flow magnitudes and the correlations between flow angles [74,99,100]. The four-particle symmetric cumulants is defined as which quantifies the correlation between v 2 m and v 2 n [101,102]. The corresponding normalized three-and four-particle cumulants are defined as nac n {3} = ac n {3}/( v 2 n · v 2 2n ), nsc m,n {4} = sc m,n {4}/( v 2 m · v 2 n ), which try to eliminate the dependence on the flow coefficients and focus on evaluating the relative strength of the correlations between different flow harmonics. Since the related calculations are numerically expansive, we only show the results nac n {3}, nsc 2,3 {4} and nsc 2,4 {4} before and after hydrodynamic evolution with TRENTo initial condition. Fig. 6 shows that nac n {3}, nsc 2,3 {4} and nsc 2,4 {4} in the final states keep the same sign of those in the initial state correlations. Another interesting feature is that the hierarchy of the four-particle correlations in final states does not follow the one in the initial states. For example, Fig. 6 (b) shows that the nsc ε 2,3 {4} 2 from three sets of parameters follows Para-I > Para-II > Para-III, but the hierarchy of the nsc v 2,3 {4} is inverted after the hydrodynamic evolutions with Para-III > Para-II > Para-I. This can be caused by the different non-linear response effects to various initial conditions. Such non-linear response effects are the greatest for Para-III, which lead to the largest nsc v 2,3 {4} after the hydrodynamic evolution. This is also consistent with 2 Here, the nsc ε 2,3 {4} and nsc ε 2,4 {4} are calculated in the 0-0.1% centrality bin in the initial states.
the results of c v 2 {4} in Fig. 3, which shows that c v 2 {4} of Para-III is the most positive one due to the largest non-linear response of v 2 .
As shown in Figs. 3 and 6, the current hydrodynamic calculations with these three initial conditions have difficulties in describing the measured multi-particle single cumulants and mixed harmonic cumulants for highmultiplicity p-p collisions. Nevertheless the presented hydrodynamic calculations also confirm that the mixed harmonic multi-particle correlations are very sensitive to the details of initial conditions. If hydrodynamics works for the small p-p collision systems, the related experimental data is very useful to constrain the corresponding initial conditions.

IV. SUMMARY
In this paper, we investigated the hydrodynamic flow in high-multiplicity events of proton-proton collisions at √ s = 13 TeV, using iEBE-VISHNU hybrid model with HIJING, super-MC and TRENTo initial conditions. With properly tuned parameters, iEBE-VISHNU can roughly reproduce the measured two-particle correlations, including the integrated and differential flow for all charged and identified hadrons. However, the hydrodynamic calculations with any initial condition can not describe the negative c 2 {4} measured in experiments, which give a wrong sign. Further investigations showed that the elliptic flow v 2 does not linearly respond to the initial eccentricity ε 2 . The non-linear (cubic) response becomes important in the small systems, which plays a non-negligible role and enhances the flow fluctuations. Such contribution always leads to a positive c v 2 {4} even when the sign of c ε 2 {4} is negative in the initial conditions.
We also performed the first hydrodynamic calculations for normalized three-and four-particle azimuthal correla-tions, nac n {3}, nsc 2,3 {4} and nsc 2,4 {4} in p-p collisions at √ s = 13 TeV, and found that iEBE-VISHNU can qualitatively describe the features of nac n {3} and nsc 2,4 {4} but fail to reproduce the negative nsc 2,3 {4} as measured in experiments. At the current stage, it is still challenging to describe the measured multi-particle cumulants of single and mixed harmonics within the framework of 2+1D hydrodynamics with these three initial conditions implemented in this paper. In the near future, it is worthwhile to implement 3+1D hydrodynamics with longitudinal fluctuations and dynamical initial conditions to further investigate these flow data in p-p collisions, which could help us to evaluate whether or not tiny droplets with collective expansion have been created in p-p collisions at the LHC.