Signatures of collective flow in high multiplicity pp collisions

A blast-wave parametrization, including a full set of hadronic resonances, is used to model a small system, with total particle multiplicity comparable to the one measured in the high-multiplicity pp collisions at the LHC. Calculations are preformed for three cases: with negligible, regular and strong radial flow on the blast-wave hypersurface. We investigate the effects of flow on inclusive p_T spectra as well as on 1D and 3D femtoscopic radii for pions. Special emphasis is put on the role of pions from resonance decays. In particular we show that they magnify the flow effects present in the blast-wave stage and significantly influence the shape of the correlation functions. A specific observable, the R^E_out/R^G_side ratio is proposed as a sensitive probe of the collective effects. Model results for the high multiplicity pp collisions, for scenarios with small and large radial flow are compared.


I. INTRODUCTION
The pp program at the LHC has successfully started, and has already produced a large set of minimum-bias measurements from the experiments [1][2][3][4][5][6][7][8]. Of particular interest is the fact that the multiplicity range measured at √ s = 7 TeV overlaps with the one measured in peripheral CuCu collisions at ultra-relativistic energies of √ s N N = 62 and 200 GeV at RHIC. This allows for comparison of observables between p+p and heavy-ion collisions at the same measured final state multiplicity dN ch /dη. In heavy-ion collisions, dN ch /dη is the scaling variable for many soft physics observables, including the inclusive particle p T spectra, pion femtoscopic radii [9], strange particle yields etc. It is therefore interesting whether a high-multiplicity pp collision resembles a heavy-ion one.
This question has become even more relevant with the recent observation of the long-range η ridge in high multiplicity pp collisions by the CMS experiment [10]. There are attempts to explain the phenomenon in the Color Glass Condensate framework [11] or as an elliptic flow [12]. Both of these critically depend on the observation of radial flow (or space-momentum correlation) in such collisions. A recent study with EPOS model [13] suggests that such flow has actually been seen in the ALICE pion femtoscopy data [4]. Similar topics have been investigated in the framework of the rescattering model [14]. A comparison of pion femtoscopic radii in pp and heavy-ion collisions at lower energy √ s = 200 GeV have been carried out by the STAR experiment [15]. Scaling between radii in pp and heavy-ion collisions has been observed vs. the pair momentum, suggesting that there might be a common physics mechanism, i.e. radial flow, driving them in both cases.
The hydrodynamic scenario is well applicable to heavyion collisions at RHIC. It produces specific space- * Adam.Kisiel@cern.ch momentum correlation patterns, which are commonly referred to as flow. The system created in the heavy-ion collision expands rapidly outwards, showing a strong radial flow, which is observed in the modification of the single particle inclusive p T spectra shape. The spectrum is an observable depending only on the momenta of the particles, so its connection to space-time can only be indirectly inferred. To directly access space-time information we employ femtoscopic techniques. The fall of the "femtoscopic radii" (sometimes called "HBT radii") with particle's m T can be interpreted as the decrease of "lengths of homogeneity", a direct consequence of radial and longitudinal flow [16]. The simultaneous description of the p T spectra and the femtoscopic radii dependence on m T in heavy-ion collisions can be achieved with the "blast-wave" parametrization [17][18][19]. In this work we test what results such an approach brings when applied to small systems. We characterize the size of the system by the average charged particle multiplicity per unit of pseudorapidity that it produces. This means that simple "blast-wave" models cannot be used, as abundances of each particle species are free parameters there. In contrast the THERMINATOR model [20] has both necessary features: it can use the traditional "blast-wave" hypersurface and provides absolute particle multiplicities.
The usage of the THERMINATOR makes it possible to test the impact of the strongly decaying resonances. As it was shown in [19], their introduction increases the apparent size of the system produced in heavy-ion collisions by approximately 1 fm and leads to the development of long-range tails in the emission function of pions. These effects were important but not dominant in the AuAu system (with the average size of 6 fm), but they may dominate in a pp collision where we expect the initial system to have a size of the order of 1 fm.
The paper is organized as follows. In Section II we briefly introduce the THERMINATOR model and its blastwave variant. In Section III we describe the choice of parameters for the three simulations that were performed. In Section IV we describe the results for the inclusive transverse momentum p T spectra and the 1dimensional femtoscopic simulations and explain the features observed there. In Section V we show the results of the 3D femtoscopic analysis and propose the most robust and sensitive observables that can be used to search for collectivity in the pp collisions at large multiplicities, e.g. at the LHC. In Section VI we discuss the origin of the effects seen in previous Sections.

II. THERMINATOR MODEL AND FEMTOSCOPIC FORMALISM
In this work we use the THERMINATOR model [20], more specifically its variant that uses the blast-wave freeze-out hypersurface with the possibility to adjust the space-time correlation at the freeze-out hypersurface. For a detailed description of the model, its full theoretical basis and formalism see [19]. In short, average per-event particle abundances are calculated from the chemical model, where the temperature T and bariochemical potential µ B are model parameters, and isospin µ I and strange µ S potential are calculated from conservation laws. For each event the number of particles of each type (we use 381 particle types taken from PDG [21]) is randomly generated from a Poissonian distribution. Then each particle is put on the freeze-out hypersurface according to the generalized blast-wave prescription [19]: where a, τ , ρ max , are parameters of the model. ρ max is the transverse size of the system, τ is the proper time of emission, and a determines the transverse position (ρ) vs. time (t) correlation at the hypersurface. Of particular importance is the form of the radial flow velocity, which depends on particle emission position: where v T is a parameter which controls the amount of flow in the system. Large values of v T give small flow, small values of v T give large flow, growing semi-linearly from 0 at the center of the source to the maximum value of 1/(1+v T ) at the edge. y, p T and ϕ are the components of the momentum of the particle (rapidity, transverse momentum and azimuthal angle respectively), while α , ρ and φ are the corresponding coordinates of the emission point. As the last step, all unstable particles are propagated and decayed, in cascades if necessary, until only stable particles remain. The emission point of the particle is either its creation point on the hypersurface (we call such particles "primordial") or a decay point of its parent resonance. Final state particle rescatterings after the emission are not included. We have estimated that for the particle densities in this work, each particle would rescatter on the average 0.8 times. Most of the rescattering points would be very close to the original emission point, so the inclusion of rescattering would not change the emission shapes significantly.
Further in the work we show results of the analysis of the identical pion femtoscopic correlations (some-times called "Bose-Einstein correlations" or "HBT correlations"). All the analyzes were performed using the standard "two-particle" method: pions from the same event were combined into pairs, and assigned a weight equal to their wave-function squared, to form the signal. The pair wave function contained only the relevant symmetrization; final state coulomb and strong interactions were ignored as for pions they are small and are not important for the determination of the femtoscopic radii. Then, pions from different events were combined into pairs to form the background. The correlation function was constructed by dividing the signal by the background. It was then fitted with various functional forms, which will be described in detail in Section IV. For an extensive description of the two-particle method and the formalism see [19].

III. SIMULATION DESCRIPTION
We performed the simulations with the blast-wave model with the following parameters: transverse size of the freeze-out hypersurface ρ max = 1.5 fm, proper time of freeze-out τ 0 = 2.0 fm/c, freeze-out temperature T f = 165.6 MeV, bariochemical potential µ B = 28.5 MeV. The freeze-out hypersurface was slightly sloped in the ρ − t plane (a = −0.2). The three simulations differed only in the value of the v T parameter; the "no-flow" had v T = 10, the "regular flow": v T = 1 and the "strong flow": v T = 0.5. The resulting flow velocity profiles can be seen in Fig. 1 as dashed lines. The common radial flow present in the model results in strong space-momentum ρ − p T correlations of the produced particles, which is reflected in the the common radial flow: where averaging is done over particles. The results are shown in Fig. 1 as symbols. The plot illustrates that the "no flow" scenario has very little ρ − p T correlations, while in the "strong flow" scenario this correlation is significant. For each of the scenarios we have simulated 2M events.

IV. RESULTS
Fixing the model parameters, as described above, fixes the average absolute multiplicity of particles per event. pions kaons protons all "strong flow" 234 ± 1 323 ± 4 414 ± 8 "no flow" 157 ± 1 168 ± 2 157 ± 3 primordial "strong flow" 285 ± 2 360 ± 6 465 ± 2 "no flow" 183 ± 1 181 ± 3 176 ± 8 The resulting dN ch /dη for the "strong flow" sample was 8.3. This is about twice the average multiplicity reported by ALICE collaboration at √ s = 900 GeV and is above the average multiplicity for √ s = 7 TeV [2][3][4]. The impact that the radial flow has on particle spectra is shown in Fig. 2. Only for the simplest case of "no flow" and primordial particles only (shown in panel d)) the spectra for all particle types are exponential. We fit the distributions with an exponential and extract their slope, shown in Tab. I. Even in the "no flow" case the simple addition of strongly decaying resonances obscures the picture, as seen in panel c). Resonances tend to populate the lower m T regions more abundantly than higher m T . There are relatively more resonances decaying into pions and protons than into kaons. As a result, the slope parameter is higher for pions and protons and lower for kaons. For primordial particles in "strong flow" in panel b) we see, as expected, that the thermal "exponential" shapes are modified by radial flow: pions get a "concave" shape, while kaons and protons develop a positive curvature. All are less steep than in the "no flow" case, resulting in the larger slope parameter. Adding resonance decay products, shown in panel a), decreases the slope parameter by 50 MeV, so any attempt to extract the flow velocity and temperature from the slope parameters must fully take into account the resonance contribution. In Fig. 3 the kaon to pion and proton to pion ratios vs. p T are shown. For K/π both flow and addition of resonances decrease the ratio. This comes from the larger contribution of resonance decays to pions. In contrast, for p/π only flow decreases the ratio, as the relative contributions of resonances to both pions and protons are similar. In summary, the particle spectra retain the general characteristics associated with radial flow even when combined with resonance decays in a small system, but resonance decays do alter the quantitative estimates of temperature significantly. In addition our simple model does not include the hard processes which will contribute to all particle's spectra at larger p T and complicate the picture further. Therefore making strong conclusions about radial flow from particle spectra alone is not trivial and is model dependent.
Another observable which was extensively used in conjunction with p T spectra as a signature of collective flow is the pion femtoscopic radius of the system, more specifically its monotonic decrease with growing pair momentum k T = |p T ,1 + p T ,2 | /2. The first measurement of such a dependence at the LHC was reported in [4]. The minimum-bias measurement shows no significant dependence, however the combined multiplicity vs. k T analysis shows that the behavior of "low" and "high" multiplicity events is different. The "low" multiplicity (M) events show an increase of pion femtoscopic radius R inv at low k T and then a fall at high k T . The "high" M events on the other hand suggest a decrease of the R inv radius with k T , although it is also consistent with a flat dependence. New results from the LHC, especially from collisions at the higher energy of √ s = 7 TeV will deliver results with better precision and will show if the trend continues at higher multiplicities.
We have obtained the pion femtoscopic correlation functions using the two-particle method described in Section II. We divide the sample into six k T bins: We begin by analyzing the 1dimensional correlation functions. Following the experiments [4,7], we fit them with two functional forms, the Gaussian: and the exponential: An example of the fits is shown in Fig. 4. We see a behavior similar to the one reported by the LHC experiments: the correlation function is better described by an exponential. The Gaussian fit is worse, but it characterizes the overall width of the correlation function. In Fig. 5 we plot the results of the 1-dimensional analysis vs. k T for the three scenarios of flow. The Gaussian and properly scaled exponential fit results are consistent. In the "no flow" scenario the R G inv is growing with k T , similar to low multiplicity collisions at the LHC. On the other hand, as flow develops, the R G inv becomes flat vs. k T for "regular flow" and has a negative slope for "strong flow". In other words, from the comparison of the R G inv vs. k T trends between data and the model, we observe a development of radial flow with the increase of multiplicity in the pp collisions at the LHC. It would be interesting to see if the data at higher multiplicities at √ s = 7 TeV continue the trend and develop an even stronger k T dependence. The results are intriguing enough to attempt a more in-depth explanation of their origin. We aim to propose more strict tests, which would make the conclusions from such comparisons less model dependent.

V. 3D BERTSCH-PRATT ANALYSIS
We have performed the 3D analysis of the simulated correlation functions. We have used the standard Bertsch-Pratt decomposition of the relative momentum into the "long" (along the beam axis), the "out" (along the pair transverse momentum), and the "side" (perpendicular to the other two) directions. We have used the Longitudinally Co-Moving System (LCMS) reference frame, where the longitudinal momentum of the pair vanishes. We have analyzed two classes of pions: the "primordial" only, and "all" particles. The latter includes the "primordial" plus those from strongly decaying resonances. In Fig. 6 we show the correlation function, separately in out, side and long directions. We attempted to fit both the 1d projections and the full 3d function with combinations of Gaussian, exponential and Lorentzian shapes. These three were chosen because they have an analytic relation between correlation function and pair emission function forms. We found that (a) for primordial particles a 3D Gaussian is the best fit, (b) for all particles the out and long direction is best fitted with an exponential, while side remains Gaussian. This is true both for 1d projections (shown in the figure) as well as for a full 3d fit. Even though the fit in Fig. 6 is best out of three forms tried, it is still not perfect. Nevertheless, the width of the function is adequately estimated, which is enough for the purpose of this paper. The investigation of the details of the shape of the correlation functions we leave for future studies. Based on the study of correlation function shapes, for primordial particles a Gaussian 3d ellipsoid, with different sizes in all directions is a reasonable description of the correlation: where R G out , R G side and R G long are Gaussian radii in out, side and long directions respectively, λ is the strength of the correlation. For all particles we have used the following fit functional form: where R E out and R E long are exponential radii in out and side directions. Similarly to the 1-dimensional case, the radii obtained from Gaussian and exponential fits are not directly comparable, as they are parameters of different functional forms.
In Fig 7 we show the Gaussian femtoscopic radii obtained from the analysis of "primordial" particles. We concentrate on the transverse dynamics. For the R G radius we observe an expected behavior: the "no flow" scenario shows a flat k T dependence, the slope grows with increasing flow and is significant for "strong flow". In addition, the size gets smaller with increasing flow; the well-knows mechanism of the "lengths of homogeneity" is at work here. R G side dependence is weakly changed by flow; stronger flow creates a small slope. The change in behavior is best visible in the R G out /R G side ratio. For "no flow" it is large, at least 1.0-1.5, for "strong flow" it is falling, and smaller, below 1.0. The R G out /R G side ratio has an advantage that it is, to a limited degree, independent of scale. The absolute values of the transverse radii are directly depending on the assumed system size (the ρ max parameter), while the ratio is less dependent on it, and its qualitative behavior should not change as we change ρ max in a range which would produce multiplicities observed in pp collisions.
For primordial particles we have identified the R G out /R G side ratio as a sensitive, scale-independent, probe of the amount of collectivity in the system. However, the majority of pions are coming from resonance decays, which strongly influence spatial and temporal characteristics of the emission process. It is therefore crucial to investigate how they modify the signal present for "primordial" particles. In Fig. 8 we show the fit results for all pions. We stress that inclusion of resonances changes the emission patterns visibly, which is reflected in the change of the fit functional form. R E out for "no flow" is larger than for "primordial" and develops a k T slope. In contrast, the "strong flow" R E out grows little: R E out grows less as flow grows. That is surprising, as the "length of homogeneity" argument seems to hold for particles from resonances, even though they do not come directly from a "collective" medium. Also the negative slope of the R E out dependence on k T alone is not a good signature of collectivity as such behavior also develops for "no flow". R G side shows the behavior opposite to R E out . For "no flow" it is similar to "primordial" case, while for "strong flow" it grows. As a consequence, again in striking contrast to the behavior of R E out , R G side grows more as radial flow increases. Again the R E out /R G side ratio is the most sensitive to these changes, and for all particles it shows an even more dramatic difference between flow scenarios: the value for the ratio is up to a factor of 1.5 larger for "no flow" as compared to "strong flow". In summary, we show that a small value of R E out /R G side ratio, close to unity, is a signal of emission of all particles from a collective medium. The value is monotonically dependent on the strength of flow: the stronger the flow, the lower the R E out /R G side ratio. In addition, the inclusion of resonance decay products magnifies this flow effect.

VI. ORIGINS OF FEMTOSCOPIC RADII BEHAVIOR
The identification of the R E out /R G side ratio as a sensitive probe of collective particle emission in systems resembling pp collisions at the LHC has consequences for experiments. It allows to test a hypothesis of a creation of a collective system in such collisions, which produces large particle multiplicities. In addition, the femtoscopic pion analysis is practically a "minimum-bias" measurement and does not require advanced triggering or acceptance corrections. It can therefore be performed quickly, with the data samples that the LHC experiments already have on tape. It is therefore critical to understand why such effect arises and how general their emergence is; in other words how model dependent is the link between R E out /R G side ratio and system collectivity. The separation between two particles in the out direction is: where x = ρp T cos(φ − ϕ)/p T are the transverse emission point coordinates of the first and the second particle, projected on the pair velocity, t 1 and t 2 are the emission times of the particles and β is the pair velocity (we ne-glect the differences between particles velocities since in order to be femtoscopically correlated their relative velocity must be small). We order the particles in a pair in such a way that t 2 > t 1 , the third term in the formula reflects the fact that the first particle has already traveled some distance before the second was born. The definition of the side direction is that the pair velocity in this direction vanishes, so the formula for the separation is simpler: where y = (ρp T sin(φ−ϕ))/p T are the transverse emission point coordinates of the first and second particle, now perpendicular to the pair velocity. x coordinates are on the x axis of the top panels of Figs. 9 and 10, y are on the x axis of the bottom panels and t is on the y axis of all panels. The R E,G out and R G side femtoscopic radii have the meaning of variances of out and side distributions respectively (conventionally divided by √ 2). They are proportional to a combination of the single particle spacetime distribution moments: where σ x , σ y and σ t are the variances of the x, y and t distributions respectively and C(x, βt) is the term that accounts for a possible space-time correlation in the emission process. We note that C contributes to the size with a negative sign: it means that positive x − t correlation decreases the out radius [22,23]. In Fig. 9 the space-time characteristics of the emission process for the "no flow" scenario is shown. On panels a) and d) are primordial particles. In this case both out and side directions show similar size, which does not depend on time. Also, the mean emission points are not shifted. Panels b) and e) show the behavior for resonance decay products. There is a small correlation between mean x and time. Also, out and side spreads grow with time, but only slightly. As a consequence, for all particles shown in panels c) and f), the out emission point is weakly correlated with time and the side spread grows only slightly.
In Fig. 10 we present a similar study for the "strong flow" scenario. We concentrate on the differences between this scenario and the "no flow" one. For primordial particles, the out mean position is shifted by flow. The spread is slightly smaller than side; hence the smaller than one ratio for green points in panel a) of Fig. 7. For resonance products we see more striking differences. The out mean position is strongly positively correlated with time. The spread also grows with time, but less than for the side, where the growth is more pronounced than for the "no flow" case. As a result, for all particles, the out mean position is positively correlated with time, while the side spread grows stronger than in the "no flow" case.
We can now understand the dramatic change of the R E out /R G side ratio between "no flow" and "strong flow" scenarios. For "no flow", σ x and σ y are similar, as can be seen in panels c) and f) of Fig. 9. But R E out is further increased by σ t and the x − t correlation is too small to counter-balance it. As a result R E out /R G side is large. In short, resonance decays introduce significant spread in emission times, which acts to increase R E out . The fact that resonance decay points are distributed according to the exponential decay time multiplied by velocity also explains why the shape of the correlation in the out direction changes dramatically, while it stays Gaussian in side, where the velocity vanishes.
In the "strong flow" scenario we note that even for primordial particles, seen in panel a) and d) of Fig. 10, σ x is smaller than σ y -the effect of flow of primordial pions. The effect is further increased by resonances products. Because of the existence of common radial flow in the system, when resonances are born on the freeze-out hypersurface, they are already flying to the outside of the system. Since they decay after some time, the emission points of the daughter pions will, on average, be further apart compared to the "no flow" case. This is nicely illustrated as the strong growth of the σ y with time on panel e) of Fig. 10. The result is a significant increase of R G side , as compared to the "no flow" case. One expects the same argument to hold for σ x , and to some extent it is seen in panel b) of Fig. 10. However, the growth is not as strong as for σ y -this is because the resonance emission points are also governed by flow, so their parent σ x is smaller than σ y , just like for primordial pions on panels a) and d). But more importantly, radial flow of resonances also results in much stronger (as compared to "no flow" case) positive x − t correlation, seen most clearly in panel b). This correlation also dominates when we consider both primordial and resonance daughter pions in panel c). This will further decrease the R E out parameter, due to the C(x, βt) term in Eq. 10. As a result, the R E out /R G side ratio is small. In short, radial flow of primordial pions and resonances produces several effects which all combine together to bring the R E out /R G side ratio down, even below unity if the flow is strong enough. The stronger the flow, the lower the ratio.
We also now understand why 1-dimensional R inv radii dependence on k T changes as flow develops. R inv is approximately equal to the quadratic average of γR out , R side and R long , where the γ factor in front of R out comes from the boost (with pair velocity) from the LCMS to Pair Rest Frame, where R inv is determined. For "no flow" case R out is flat with k T in LCMS, so the growth of the boost as k T increases drives R inv up. For cases with flow, the growth of the boost is countered by the fall of R out resulting from flow. And since the other two radii also decrease with k T , the slope of the R inv dependence on k T becomes negative.
As a final note we stress that the model employed here is a simple hydrodynamics-inspired parametrization. We have shown that one of its particular features -an ρ − p T correlation, also known as radial flow, produces specific signatures in single particle spectra as well as 1D and 3D pion femtoscopic radii. However it does not address the question of where such correlations come from. One possibility is that they come from a "medium" that is described by hydrodynamic equations, as suggested in [13], but there may be other mechanisms that produce such correlations.

VII. CONCLUSIONS
We have performed simulations with the THERMINATOR model, employing emission of primordial particles from a blast-wave hypersurface, and including the propagation and decay of a full set of hadronic resonances. The parameters of the calculations have been adjusted, so that they produce multiplicities observed in high multiplicity pp collisions at LHC energies. Three sets of parameters have been used, which differ only by the strength of the radial flow (the ρ − p T correlation) in the system.
The radial flow modifies the particle spectra in the expected way, creating a concave shape for pions and positive curvature for kaons and protons. However, addition of resonance decays significantly modifies these shapes, influencing the extracted slope parameter and making a direct extraction of flow velocity from spectra measurements alone complicated and model dependent.
Results of 1-dimensional femtoscopic analysis show that with no flow, the 1D R G inv pion femtoscopic radius is expected to grow with pair momentum k T . In contrast, as flow is increased, the dependence changes to a decrease with k T . Therefore, within the frame of this model the results of ALICE collaboration on multiplicity vs. k T dependence can be interpreted as signature of the development of radial flow with increasing per-event particle multiplicity.
3D pion femtoscopic radii show that the shape of the correlation is exponential in the out and long directions due to contributions from resonance decay products, while in side the shape remains Gaussian. Radial flow (or lack of it) for primordial pions and resonances has opposite effect on R E out and R G side radii. The first grows less as flow develops, which is a result of a combination of two effects: (a) decrease of lengths of homogeneity for primordial pions and resonances and (b) strong x − t correlations coming from resonances' flow-ordered velocities. R G side grows more as flow develops, the effect of the preferential propagation of the resonances to the outside of the source in the presence of flow. We propose to use the R E out /R G side ratio as a probe of the amount of radial flow (e.g. velocities ordering) in the small, resonancedominated system, produced in the high multiplicity collisions at the LHC. Any analysis aiming to search for collective effects in such systems should combine at least the three measurements: identified particle p T spectra, 1D and 3D femtoscopic radii.