1 Introduction

Ultra-relativistic collisions of heavy ions are intended to create a novel state of matter, the Quark-Gluon Plasma (QGP), and to study its properties. Extensive measurements of various flow observables performed at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) together with the successful descriptions from hydrodynamic calculations revealed that the created QGP fireball behaves like a nearly perfect liquid with very small specific shear viscosity [1,2,3,4,5,6,7,8,9,10,11,12,13,14]. Recently, various striking features of collective expansion have been observed in high-multiplicity events of the small collision systems, such as p–Au, d–Au, \(^3\)He–Au at RHIC [15, 16] and p–p and p–Pb at the LHC [17,18,19]. These features include the long-range “double ridge” structures in two-particle azimuthal correlations with a large pseudo-rapidity gap even up to 8 units [20,21,22,23,24,25,26,27,28,29,30,31,32], the changing signs of the 4-particle cumulants [15, 27,28,29,30,31,32,33,34] and \(v_2\) mass ordering of identified hadrons [30, 35,36,37], etc.

These observed flow-like signals in the small systems can be quantitatively or semi-quantitatively described by 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,57,58,59,60,61], hadronic rescatterings [62, 63], a string rope and shoving mechanism [64] 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 [65,66,67,68,69,70,71,72,73,74] 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 [75] 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 \(\sqrt{s_{NN}}=5.02\) TeV, which demonstrate the importance of the partonic degrees of freedom and possible formation of QGP in the small p–Pb systems.

Recently, the collectivity and possible formation of QGP in high-multiplicity p–p collisions at the LHC energies has also attracted lots of attention. Compared to larger collision systems, the corresponding non-flow contributions, such as mini-jets or resonance decays, become more significant. In the measurements of two-particle correlations, two different non-flow subtraction methods, template fit [24, 29, 31] and peripheral subtraction [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 short-range two-particle correlations [76, 77]. Besides, two- and 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 [77,78,79]. It was found that \(c_2\{4\}\) turns to negative value in high-multiplicity 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 [78, 79]. Furthermore, ALICE [32], ATLAS [79] and CMS [80] 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 [81], we found that, with properly tuned parameters, hydrodynamic simulations with HIJING initial conditions can nicely describe the two-particle correlations in p–p collisions at \(\sqrt{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 \(\varLambda \)). 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 [82] and TRENTo [83], 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 four-particle azimuthal correlations. To better understand the non-linear hydrodynamic evolution in the small systems, we also investigate the response between the initial \(\varepsilon _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 hydrodynamic model and the initial conditions of HIJING, super-MC and TRENTo and explain the setups. Section 3 presents the model calculations, the comparison to the experimental data, and the related discussion. Section 4 gives a brief summary of this paper.

2 The model and set-ups

2.1 iEBE-VISHNU hybrid model

iEBE-VISHNU [84] is an event-by-event version of hybrid model VISHNU [85] that combines 2+1D viscous hydrodynamics VISH2+1 [86, 87] to describe the QGP expansion with a hadron cascades model UrQMD [88, 89] 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^{\mu \nu }\) and shear stress tensor \(\pi ^{\mu \nu }\) with a state-of-the-art equation of state (EoS) s95-PCE [90, 91] 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 \(\eta /s\). The hydrodynamic evolution matches the hadron cascade simulations at a switching temperature \(T_{\mathrm{sw}}\), where various hadrons are emitted from the switching hyper-surface 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 [82] and TRENTo [83]. 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.

2.2 HIJING initial condition

In HIJING  [92,93,94], 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 \(\sigma _R\): \(\exp [-\frac{(x_i-x_{c})^{2}+(y_i-y_{c})^{2}}{2\sigma _{R}^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]

$$\begin{aligned} e(x,y)= & {} K\sum _{i}\frac{p_{i}U_{0}}{2\pi \sigma _0^{2}\tau _{0}\varDelta \eta _{s}}\nonumber \\&\exp \left[ -\frac{(x-x_{i})^{2}+(y-y_{i})^{2}}{2\sigma _0^2}\right] , \end{aligned}$$
(1)

where \(\sigma _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 \(|\eta |<1\) (for related details, please also refer to Ref. [51]).

Table 1 Four parameter sets of iEBE-VISHNU simulations with HIJING initial condition for p–p collisions at \(\sqrt{s}=\) 13 TeV
Full size table

2.3 super-MC initial condition

For p–p collisions, super-MC model with sub-nucleonic fluctuations [82] 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. [82], the initial entropy density of the produced matter is modeled as

$$\begin{aligned} s({\mathbf {r}}) = \frac{\kappa _s}{\tau _0} \sum _{k=1}^3 \frac{\gamma _k^{(i)}}{2\pi \sigma _g^2} \exp \biggl [-\frac{({\mathbf {r}}-{\mathbf {r}}_k^{(i)})^2}{2\sigma _g^2}\biggr ], \end{aligned}$$
(2)

where \(\gamma _k^{(i)}\) (\(i=1,2,3\)) are the random weighting factors that used to fit the multiplicity distributions in p–p collisions, \(\mathbf{{r}}_i\) (\(i=1,2,3\)) is the position of three valence quarks which is distributed according to a Gaussian probability distribution. \(\sigma _g\) is a factor to describe the shape of quark density distribution together with a consideration of low-x gluon contributions [82].

2.4 TRENTo initial condition

TRENTo is a parameterized initial condition model, which generates the initial entropy density via the reduced thickness function [95, 96]:

$$\begin{aligned} s = s_0 \left( \frac{{\tilde{T}}_A^p + {\tilde{T}}_B^p}{2} \right) ^{1/p}, \end{aligned}$$
(3)

where \({\tilde{T}}(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. [83, 95, 96].

For proton–proton collisions, TRENTo is modified with the sub-nucleonic structure [83] so that \({\tilde{T}}(x, y)\) is written as \({\tilde{T}}(x, y)\equiv \int dz\, \frac{1}{n_c} \sum _{i=1}^{n_c} \gamma _i\, \rho _c \,({\mathbf {x}} - {\mathbf {x}}_i \pm {\mathbf {b}}/2)\), where \(n_c\) is the number of independent constituents in a proton, \( \gamma _i\) (\(i=1,2, ..., n_c\)) is a random weighting factor with the unit mean and variance 1/k, \({\mathbf {x}}_i\) (\(i=1,2, ..., n_c\)) are the positions of constituents, \(\mathbf {b}\) is the impact parameter, and \(\rho _c\) is the density of constituents written in a Gaussian form: \( \rho _c({\mathbf {x}}) = \frac{1}{(2\pi v^2)^{3/2}} \exp (-\frac{{\mathbf {x}}^2}{2 v^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. [83, 97, 98], 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 \(\tau _0\) to obtain the boost-invariant energy-momentum tensor \(T^{\mu \nu }(x,y, \tau _0)\). After that, we implement the following Landau matching condition to obtain the initial energy density \(e(x,y, \tau _0)\) and fluid velocity \(u^\mu (x,y, \tau _0)\):

$$\begin{aligned} T^{\mu \nu } u_\nu = e u^\mu , \end{aligned}$$
(4)

and the initial shear stress tensor and bulk pressure can be calculated with:

$$\begin{aligned}&\varPi = -\frac{1}{3} \varDelta _{\mu \nu } T^{\mu \nu } - P, \end{aligned}$$
(5)
$$\begin{aligned}&\pi ^{\mu \nu } = T^{\mu \nu } - e u^\mu u^\nu + (P + \varPi ) \varDelta ^{\mu \nu }, \end{aligned}$$
(6)

with the spatial projector being \(\varDelta ^{\mu \nu } = g^{\mu \nu } - u^\mu u^\nu \), together with an equation of state of \(P=\frac{1}{3} e\) for the massless ideal gas at the initial state.

Table 2 Three parameter sets of iEBE-VISHNU simulations with super-MC initial condition for p–p collisions at \(\sqrt{s}=\) 13 TeV
Full size table
Table 3 Three parameter sets of iEBE-VISHNU simulations with TRENTo initial condition for p–p collisions at \(\sqrt{s}=\) 13 TeV
Full size table

For p–p collisions at \(\sqrt{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 [99] and \(v_2\{2\}\) [30, 31, 79, 80] 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 1, 2 and 3. For TRENTo initial condition with the parameter set Para-I, we consider two cases with and without free-streaming as described above.

Fig. 1
figure 1

\(v_2\{2\}\), \(v_3\{2\}\) and \(v_4\{2\}\) as a function of \(N_{\mathrm{ch}}\) in p–p collisions at \(\sqrt{s}=\) 13 TeV, calculated by iEBE-VISHNU with HIJING (a), super-MC (b) and TRENTo (c) initial conditions. The CMS and ATLAS data are taken from Refs. [30, 80] and Refs. [31, 79], respectively

Full size image
Fig. 2
figure 2

\(v_2(p_T)\) for all charged hadrons (a–c), for \(K_S^0\) and \(\varLambda \) (d–f) in p–p collisions at \(\sqrt{s}= 13\) TeV, calculated by iEBE-VISHNU with HIJING, super-MC and TRENTo initial conditions. The CMS and ATLAS data are taken from Refs. [29, 30], respectively

Full size image
Fig. 3
figure 3

\(c_2\{4\}\) as a function of \(N_{\mathrm{ch}}\) in p–p collisions at \(\sqrt{s}= 13\) TeV, calculated by iEBE-VISHNU with HIJING (a), super-MC (b) 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, 31], respectively

Full size image
Fig. 4
figure 4

Event-by-event \(\varepsilon _2\) distributions \(P(\varepsilon _2)\) of HIJING (a), super-MC (b) and TRENTo (c) initial conditions at 0–0.1% centrality bin in p–p collisions at \(\sqrt{s}=13\ \text {TeV}\)

Full size image
Fig. 5
figure 5

Left panel: the scatter points between the \(v_2\) and \(\varepsilon _2\), together with a linear fitting and a non-linear fitting with both linear and cubic terms. Right panel: the comparison between the scaled event-by-event \(\varepsilon _2\) distribution and scaled \(v_2\) distributions for iEBE-VISHNU simulations with HIJING initial condition (Para-III) at 0–0.1% p–p at \(\sqrt{s}= 13\) TeV

Full size image

3 Results and discussion

3.1 2-Particle cumulant

With various sets of parameters for these three initial conditions, HIJING, super-MC and TRENTo, as listed in Tables 1, 2 and 3, we calculate the integrated \(v_n\{2\}\) (\(n=2\), 3 and 4) as a function of multiplicity for p–p collisions at \(\sqrt{s}=\) 13 TeV, using iEBE-VISHNU together with an application of the two-subevent method with the pseudorapidity gap \(|\varDelta \eta | > 0\), kinematic cuts \(0.3< p_{\mathrm{T}} < 3.0 \) GeV/c and \(|\eta |<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_{\mathrm{ch}}^{\mathrm{sel}}\) with \(0.3< p_{\mathrm{T}} < 3.0 \) GeV/c and \(|\eta |<2.4\), and then map it to the number of charged hadrons \(N_{\mathrm{ch}}\) with \(0.4 < p_{\mathrm{T}} \) GeV/c and \(|\eta |<2.4\) to compare with the experimental data. Figure 1 presents the comparison between our hydrodynamic calculations and the experimental measurements from ATLAS [30, 80] and CMS [31, 79]. 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 [79].

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), which free-streams the initial state to proper time \(\tau _0\) before instantaneously switching to hydrodynamic simulations [83, 98]. 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. Such difference of the calculations with and without free streaming at the pre-equilibrium stage is mainly caused by the treatment of the initial flow and the initial shear-stress tensor. In the model calculations without the pre-equilibrium dynamics, we set the initial flow and shear stress to zero when the hydrodynamic start. While the calculations with pre-equilibrium dynamics include the initial flow and shear stress calculated by the energy-momentum tensor of the particles, the qualitative behavior in the right panel of Fig. 1 (red dashed line) changes in the calculations with the pre-equilibrium dynamics. Its slope of the \(v_2\{2\}\) as a function of \(N_{ch}\) becomes larger. This can be understood in the following way:

The higher multiplicity events will have higher initial energies density as well as the larger gradients of the energies. Thus, these events have stronger initial flows, which enlarges the final \(v_2\{2\}\) compared to lower multiplicity events. Note that the lifetime of hydrodynamic evolution is short in p–p system, so the initial flow will have relatively larger effects on the final observables compare to the case of the large collision systems. For the more studies about the effects of the pre-equilibrium dynamic evolution described by the kinetic model, please refer to papers [52,53,54,55,56,57, 61].

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_S^0\) and \(\varLambda \) (d)–(f) for the multiplicity range \(80< N_{\mathrm{ch}}^{\mathrm{sel}} <120\) with the two-particle cumulant method with a pseudorapidity gap \(|\varDelta \eta |>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“template 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_S^0\) and \(\varLambda \) 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_s^0\) and \(\varLambda \), and the calculations with TRENTo initial condition tend to underestimates the data of \(\varLambda \). The mass splitting between \(K_S^0\) and \(\varLambda \) is more significant for calculations with TRENTo initial condition. Such larger mass splitting of \(v_2\) indicates a stronger radial flow development 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 [100]), which can also provide certain constrains on the initial conditions.

3.2 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 \(\sqrt{s} = 13\ \text {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^\varepsilon _2\{4\}\) for 0–0.1% events in the initial states, as listed in Table 4. 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 \(\sqrt{s}=13\ \text {TeV}\) [101]. We thus emphasize that hydrodynamic simulations do not necessarily produce negative \(c^v_2\{4\}\), and the observed negative \(c^v_2\{4\}\) in experiments does not necessarily suggest hydrodynamic flow in small systems.

Table 4 \(c^{\varepsilon }_2\{4\}\) for 0–0.1% centrality calculated by HIJING, super-MC, and TRENTo initial conditions with three or four sets of parameters
Full size table

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 \(\varepsilon _2\), the scaled \(v_2\) distributions \(P(v_2/\langle v_2\rangle )\) should overlap with the scaled \(\varepsilon _2\) distributions \(P(\varepsilon _2/\langle \varepsilon _2\rangle )\), which is the case for central and mid-central Pb–Pb collisions [102]. 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^{\varepsilon }\{4\}\). However, hydrodynamic simulations did not confirm such expectation. As shown in Figs. 3 and 4, and Table 4, the negative initial \(c^{\varepsilon }_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^{\varepsilon }_2\{4\}\) of Para-III of TRENTo initial conditions is more negative than that of Para-I in Table 4, the initial \(\varepsilon _2\) distribution of Para-III of TRENTo initial condition is “wider” with larger mean value of \(\varepsilon _2\). The corresponding larger non-linear effects during the evolution lead to a larger positive value of \(c_2^v\{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 \(\varepsilon _2\) and final elliptic flow \(v_2\). As shown in Fig. 5a, a clear deviation of elliptic flow from linear scaling is observed for \(\varepsilon _2 > 0.5\) where the cubic term becomes significant, which is similar to the peripheral Pb–Pb collisions [103] Footnote 1. Such non-negligible cubic response leads to the fact that the scaled distribution \(P(v_2/\langle v_2\rangle )\) and \(P(\varepsilon _2/\langle \varepsilon _2\rangle )\) does not overlap with each other as shown in Fig. 5b. It also introduces additional fluctuations of \(v_2\) in the final states, which could even change the sign of \(c_2^v\{4\}\) and make the model calculations fail to reproduce the negative \(c_2^v\{4\}\) measured in experiments.

Fig. 6
figure 6

\(nac_{2}\{3\}\), \(nsc_{2,3}\{4\}\) and \(nsc_{2,4}\{4\}\) as a function of \(N_{\mathrm{ch}}\) in p–p collisions at \(\sqrt{s}=\) 13 TeV, calculated by iEBE-VISHNU with TRENTo initial conditions, using standard cumulant method. \(nsc^\varepsilon _{2,3}\{4\}\) and \(nsc^\varepsilon _{2,4}\{4\}\) of the initial state in 0–0.1% centrality bin are also shown. The ATLAS data with three-subevent method are taken from Ref. [79]

Full size image

It has been generally argued that two- and multi-particle cumulants have different sensitivities to the flow fluctuations, which is written as [76]

$$\begin{aligned} v_{n}\{2\}^2= & {} \langle v_n \rangle ^2 + \sigma _v^{2}, \nonumber \\ v_{n}\{4\}^2= & {} \langle v_n \rangle ^2 - \sigma _v^{2}. \end{aligned}$$
(7)

Here \(\langle v_n \rangle \) and \(\sigma _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 \(\langle v_n \rangle \) and the flow fluctuations \(\sigma _v\).

On the other hand, the expectation of the “linear response” relation between the initial \(\varepsilon _2\) and the final \(v_2\) is purely an empirical one in the realistic setups. The linear-response relation can be broken by various effects. For example, the value of \(\varepsilon _n\) is calculated using the entire region of the initial conditions while the hydrodynamics are only solved in the region above the switching temperature. In the small systems, the hydrodynamic region becomes relatively smaller, so it is non-trivial to expect the relation. Besides, the formation of the initial \(\varepsilon _2\) might be ill defined when \(\varepsilon _2\) is large. Therefore, to describe the observed negative \(c_2\{4\}\) in data, it’s necessary to do the systematic studies for various possible mechanisms using more sophisticated models, such as the dynamical initialization model [?] combined with the kinetic equations including the non-hydrodynamic modes [54, 57].

In Fig. 6, we further study the normalized three- and four-particle azimuthal correlations in high-multiplicity proton–proton collisions at \(\sqrt{s}=\) 13 TeV. The three-particle asymmetric cumulant is defined as \(ac_n\{3\} = \langle v_{n}^2 v_{2n} \cos 2n(\varPsi _n-\varPsi _{2n})\rangle \), which is sensitive to the correlations between flow magnitudes and the correlations between flow angles [79, 104, 105]. The four-particle symmetric cumulants is defined as \(sc_{m,n}\{4\} = \langle v_{m}^2v_n^2 \rangle -\langle v_{m}^2 \rangle \langle v_{n}^2 \rangle \), which quantifies the correlation between \(v_{m}^{2}\) and \(v_{n}^{2}\) [106, 107]. The corresponding normalized three- and four- particle cumulants are defined as \(nac_n\{3\}=ac_n\{3\}/(\langle v_{n}^2\rangle \cdot \sqrt{\langle v_{2n}^2 \rangle })\), \(nsc_{m,n}\{4\} = sc_{m,n}\{4\}/(\langle v_{m}^2\rangle \cdot \langle v_n^2 \rangle )\), 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. Figure 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. 6b shows that the \(nsc_{2,3}^\varepsilon \{4\}\) Footnote 2 from three sets of parameters follows Para-I > Para-II > Para-III, but the hierarchy of the \(nsc_{2,3}^v\{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_{2,3}^v\{4\}\) after the hydrodynamic evolution. This is also consistent with the results of \(c_2^v\{4\}\) in Fig. 3, which shows that \(c_2^v\{4\}\) of Para-III is the most positive one due to the largest non-linear response of \(v_2\).

Last but not least, it should be emphasised that hydrodynamic calculations, no matter which set of parameters are used, can not reproduced the measurements shown in Fig. 6, it is even worse in the case of \(nsc_{2,3}^v\{4\}\) where the sign of hydrodynamic calculation is opposite with data. 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 high-multiplicity 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.

4 Summary

In this paper, we investigated the hydrodynamic flow in high-multiplicity events of proton–proton collisions at \(\sqrt{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 \(\varepsilon _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^\varepsilon _2\{4\}\) is negative in the initial conditions.

We also performed the first hydrodynamic calculations for normalized three- and four-particle azimuthal correlations, \(nac_n\{3\}\), \(nsc_{2,3}\{4\}\) and \(nsc_{2,4}\{4\}\) in p–p collisions at \(\sqrt{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. It is found that the critical challenge in describing experimental measurements on multi-particle cumulants, also persist in the cases of mixed harmonics like nsc and nacs. 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. It would also be important to investigate further the effect of pre-equilibrium dynamics by quantifying the effect of hydrodynamic and non-hydrodynamic modes in p–p collisions as suggested in Ref. [54].