Study of collective anisotropies $v_2$ and $v_3$ and their fluctuations in $pA$ collisions at LHC within a relativistic transport approach
Yifeng Sun Salvatore Plumari Vincenzo Greco
nucl-th






sunyfphy@lns.infn.it
Laboratori Nazionali del Sud, INFN-LNS, Via S. Sofia 62, I-95123 Catania, Italy


salvatore.plumari@ct.infn.it
Department of Physics and Astronomy, University of Catania, Via S. Sofia 64, 1-95125 Catania, Italy
Laboratori Nazionali del Sud, INFN-LNS, Via S. Sofia 62, I-95123 Catania, Italy


greco@lns.infn.it
Laboratori Nazionali del Sud, INFN-LNS, Via S. Sofia 62, I-95123 Catania, Italy
Department of Physics and Astronomy, University of Catania, Via S. Sofia 64, 1-95125 Catania, Italy



             


We have developed a relativistic transport approach at fixed η/s(T) that incorporates initial space fluctuations generated by wounded quark model to study the hadron observables in 5.02 TeV p+Pb collisions. We find that our approach is able to correctly predict quite well several existing experimental measurements assuming a matter with η/s=1/4π, a result similar to previous studies within a viscous hydrodynamics approach.  Besides, we further discuss the sensitivity of the results on both η/s(T) and the smearing width. Our transport approach has the possibility to include in initial conditions the power law tail associated to minijet, and this improvement extends the agreement with the experimental data to higher p_T ranges. We also perform a comparison to Pb+Pb collisions pointing out that even if the collective flows have a similar magnitude the features of the matter created are different. By studying the correlation between collective flows and initial geometry, we find that the correlation decreases faster in small systems with the increase of n and centrality. In particular we show that the variance of σ_v_n/⟨ v_n⟩ has a quite different evolution with centrality for p+Pb, so their measurement could provide some further hint about the correctness of current modelling.









 



t







                              


Study of collective anisotropies v_2 and v_3 and their fluctuations in pA collisions at LHC within a relativistic transport approach
    Vincenzo Greco
    December 30, 2023
====================================================================================================================================





§ INTRODUCTION

High energy nucleus-nucleus collision experiments at the BNL Relativistic Heavy Ion Collider (RHIC)
 <cit.> and the CERN Large Hadron Collider (LHC) <cit.> have provided the convincing evidence for the formation of a quark-gluon plasma (QGP) in the early stages of a heavy ion collision. This QGP has been found to be strongly correlated and to exhibit strong collective behavior. The theoretical calculations within viscous hydrodynamics <cit.> and transport approach <cit.> have shown that such behavior is consistent with a matter with low shear viscosity to specific entropy ratio η/s ≃  0.1-0.2 close to the conjectured lower bound η/s=1/4π for a strongly interacting system <cit.>.

In recent years, experimental measurements from collisions of protons and deuterons with heavy nuclei reveal that high multiplicity events also show strikingly similar collective behavior as those seen in heavy ion collisions <cit.>. Several theoretical studies based on hydrodynamics <cit.> and transport approach <cit.> show that this behavior can also be attributed to final state interactions (hydro-like expansion)  similarly to nucleus-nucleus collisions. 

We discuss here the results obtained within a relativistic transport approach where the collision integral is tuned to 
describe the evolution for a fluid at finite η/s. Such an approach has been shown to reproduce the viscous hydrodynamics
system expansion in AA collisions <cit.> and even the ideal hydrodynamics behavior  in the infinite cross section limit <cit.>. The advantage of the transport approach is that one can naturally include initial condition with a power law tail at p_T> 2 GeV/c representing the minijet distribution and follow the dynamics also at increasing 
Knudsen number (or large η/s). In this paper we in particular show a first study, to our knowledge, of the 
the correlations between the final collective flows and initial geometry, which can shed light on the understanding if the experimental measurements are due to final state interactions or some initial condition effects <cit.>. For this end, we thus develop an event-by-event transport approach that incorporates initial fluctuations to study the collective behaviors in p+Pb collisions, and the correlations between initial geometry and final collective flows. We have also found
in particular that the variance of v_2 and v_3 as a function of centrality should show a pattern
quite different with respect to the one in Pb+Pb collisions.

This paper is organized as follows. In the next section, we will detail the model setup of our kinetic approach, and how we  generate initial conditions event-by-event. In the section following, we will compare our results on hadron observables to experimental results after constraining the parameters related to the initial conditions. Sec. IV will discuss the correlation coefficient that characterizes the correlations between collective flows and initial geometry and the distribution of them. Finally, conclusions and discussions will appear in Sec. V. 



§ MODEL SETUP


To model the evolution of the small collision systems, we are employing the same relativistic transport code developed in these years <cit.>, which has been used to study the dynamics of heavy-ion collisions at both RHIC and LHC energies, by solving the Relativistic Boltzmann Transport (RBT) equation. Numerically we solve the RBT equation using the test particle method, and the collision integral is solved by Monte Carlo method based on stochastic interpretation of transition amplitude <cit.>.

    p_μ ∂^μ f(x,p)+ m^* ∂^μ m^* ∂_p^μf(x,p)=  C [f]

The key aspect of our approach is to gauge the collision integral C [f(x,p)] locally to a specific value of 
the viscosity to  entropy density ratio η/s  <cit.>. Furthermore we can switch smoothly it to an increasing
η/s one, reaching the estimated value in hadronic phase <cit.>, when the system reaches the cross over region. 
This is realized by determining the total isotropic cross section according to the Chapmann-Enskog approximation:

    η = f(z) T/σ_tot

where z=m/T while the function f(z) is defined by the following expression

    f(z) = 15/16[ z^2 K_3(z)]^2/(15z^2+2)K_2(2z)+(3z^3+49z)K_3(2z)

where K_n-s are the modified Bessel functions. 
The entropy density for a massive system is given by s=ρ( 4 + z K_1(z)/K_2(z) ).

In the ultra-relativistic limit z→ 0 and the function f(z) → 1.2 and we recover the massless limit for the η. 
As shown in <cit.> 
The expression for η in Eq.(<ref>) is in quite good agreement
at level of 3-5% with the Green-Kubo formula <cit.>. The final hadron is recovered at the end of the evolution using parton-hadron duality ansatz.

For initial conditions of partons, we use a modified Monte Carlo Glauber model assuming three constituent quarks localized within each nucleon inspired by wounded quark model <cit.>, which can naturally obtain the linearity between the multiplicity of charged hadrons and the number of wounded quarks. 

For that reason, we firstly randomly place the constituent quarks in the nucleons, where the positions of nucleons in Pb nuclei are distributed according to the standard Woods-Saxon distribution with parameters R_0=6.5 fm and a=0.54 fm, according to the distribution dN/dr=r^2/r_0^3e^-r/r_0 with r_0=0.3 fm <cit.>, and then shift the center of mass of the three quarks to the position of the nucleon. After that, we generate the wounded quark profile using Monte Carlo Glauber model, where we decide whether each quark pair from target and projectile can collide or not with a probability p=e^-π r^2/σ_qq with σ_qq=13.6 mb in 5.02 TeV p+Pb collision <cit.>.

For the distribution of the spatial rapidity, we take the profile from Ref. <cit.>

    ρ_L±(η)=(1±η/η_m)exp(-(|η|-η_0)^2/2σ_η^2θ(|η|-η_0)),

where the parameters are chosen as η_m=5.7, η_0=2.5 and σ_η=2.5 to have the same shape of dN_ch/dη measured by ATLAS Collaboration <cit.> as shown in Fig. <ref>. This profile can account for more particles produced in the direction of nucleus, and greater asymmetry in events with larger multiplicity. 



After using the above procedure, the total initial parton density is then given as 

    dN/d^2𝐱_⊥dη=∑_i=1^N_partn_i ρ_⊥(𝐱_⊥-𝐱_𝐢)ρ_L±(η),

where 𝐱_𝐢 is the transverse position of each participant quark, n_i is the number of partons generated by each participant, and ρ_⊥(𝐱_⊥)=1/2πσ^2e^-𝐱_⊥^2/2σ^2. For the Gaussian distribution of the parton transverse density, we will change the parameter σ in the range of 0.4-0.6 fm in order to 
study its effect on the integrated v_2 at each centrality. 

As the number of particles produced in p+p collisions fluctuates according to a negative binomial distribution (NBD), we thus take n_i in Eq. (<ref>) to be n_0N, where N is sampled according to NBD P(N)=Γ(N+κ)n^Nκ^κ/Γ(κ)N!(n+κ)^N+κ, and n_0 is a constant such that the final charged particle multiplicity is same as that measured in experiments. We find that κ=0.54, n=3.9 and n_0≈2.352 can almost reproduce the distribution of charged particle measured by CMS Collaboration <cit.>, as shown in Fig. <ref>. We also mention that our N_ch is assumed to be equal to N_trk^offline after the efficiency corrections usually adopted in hydro calculation <cit.>. This however may introduce some uncertainty in selecting the beam centrality when comparing to data.
We notice that the values chosen are the same as the ones employed in an early approach to pA collisions <cit.>.



To get the momentum distribution of initial partons, we employ a blast wave model without initial transverse flow for it:

    dN/d^2𝐱_⊥dη=gτ_0/(2π)^2p_Tm_T
          ×cosh(η-y)e^-m_Tcosh(η-y)/T(𝐱_⊥,η)dp_Tdy,

where g=2×8+3×2×6=52 is the degree of freedom of partons (three flavor quarks and gluons), τ_0=0.4 fm/c is the thermalization time (again taken as in standard hydro approach for pA), and m_T=√(m^2+p_T^2) is the transverse invariant mass with m=0.3 GeV. 
This value is chosen because it entails a correct asymmetry in the pseudorapidity charged particle distribution, as shown in 
Fig. <ref>. Also in thermal equilibrium it generates an equation of state close to the one calculated in lattice QCD <cit.>.
We use the above relation to calculate the temperature locally from initial parton density, and then sample the momenta of partons according to Eq. (<ref>).

Finally, in order to compare our results to experimental data, we need to shift the rapidity of all charged hadrons from y^* in center of mass frame to y=y^*-0.465 in lab frame, if we define the movement of Pb nuclei as the positive z direction, as the beam energies were 4 TeV for protons and 1.58 TeV per nucleon for lead nuclei in 5.02 TeV p+Pb collisions.



§ COMPARISON TO EXPERIMENT



Before comparing our results directly to measured data, we firstly need to check the convergence of our simulation. To do this, we use the modified Glauber model to generate the same initial conditions of partons, and then evolve the systems until kinetic freeze out.
This is the first time we employ our RBT approach to study p collisions.
In fact the dimension and the densities explored by a pA system are  significantly smaller than the ones in AA,
therefore we performed a convergency test to fix the grid size and the N_ test appropriate to compute observables, where N_ test is the number of test particle for each real particle,
and in particular v_2 and v_3, in pA collisions.


We have found that  N_test=4000 and Δ x=Δ y=0.15 fm in the simulations guarantee a convergency 
for v_2,3(p_T) up to p_T ∼5 GeV/c.


In Fig. <ref>, we compare our spectra to experimental measurements in minimum-bias 5.02 TeV p+Pb collisions from ALICE Collaboration <cit.>, where the spectra is calculated in the center of mass frame of p+Pb collisions in both our studies and experiments. It is seen that our calculations agree with experimental ones in the range 0.3<p_T<1.5 GeV/c, while they underestimate the charged particle multiplicity at both low and high p_T due to two different reasons. The underestimation of charged particles at low p_T can be attributed to the missing of resonance decays in our approach, which is absent in this study using simple parton-hadron duality ansatz. At higher p_T  we do not get a good description of the spectrum, however this feature is 
similar to the one in hydro approaches <cit.>. On the other hand the transport approach can be naturally extended including an initial
non-equilibrium distribution with the power law tail at increasing p_T associated to the production of minijets.
We will see in the next paragraph that indeed the inclusion of minijets will allow to extend the validity region of the present approach
with respect to hydrodynamics, see Fig.<ref>.






The significantly large anisotropies v_2,3 observed experimentally in p+Pb collisions had initially be seen as surprising,
because in AA collision they were associated to formation time of about 4-5 fm/c while the high density state in p+Pb
collision is expected to be quite shorter. However in <cit.> it was already discussed 
that the v_2 in an ideal hydrodynamical
expansion (zero viscosity) with a conformal equation of state can be expected to scale with the size of the formed system at different centralities.
We therefore investigated how the v_2 and v_3 evolve with time scaling the last with the mean square radius 
R=√(⟨ r^2⟩), where ⟨ r^2⟩ is the mean square radius of the system weighted with the local density of the system.
Indeed the value of R goes from values of about 1 fm for p+Pb (almost independent on centrality) up
 to value of about 3.4 fm for Pb+Pb at  20-30% and 4.6 fm in Pb+Pb at 0-1% centrality.
In Fig. <ref> we show by red lines the time evolution of the building up of v_2 
(normalized to its asymptotic maximum value) in the upper panel in central collisions and in the lower panel in mid peripheral collisions for p+Pb (solid line) and
Pb+Pb (dashed line) and for v_3 (blue line).
We see that there is an approximate scaling with the time normalized to the size of the system R, and within
a time of about 1.5 R for v_2 and v_3 the anisotropies are already completely developed.


Using the width σ=0.55 fm for gaussian fluctuations of the parton transverse density and η/s=1/4π, we find also a good agreement with experimental data measured by CMS Collaboration <cit.> for elliptic and triangular flows, as shown in Fig. <ref>. It is seen by the solid red and blue lines that v_2 and v_3 agrees with experimental data up to p_T ≃ 2.5 GeV/c in our calculations which is quite similar to the results obtained within viscous hydrodynamics with very similar 
initial conditions <cit.>. We notice that we will compare our results with v_2{2} after jet contribution subtracted. We recall that experimentally v_2{2} is significantly larger than v_2{4}. Such an aspect should be considered more carefully in future studies especially in view of assessing to what extent the expansion of the matter created in pA is of pure hydrodynamical nature as well as for a more precise determination of η/s. The disagreement between our calculations and experimental measurement at higher p_T can be attributed to the missing of minijet production in higher transverse momentum, where non-equilibrium effect becomes more important. We also discuss this in next section. 



We also study the effect of η/s on v_2 and v_3, which is shown by the dashed red and blue lines in Fig. <ref>. It is seen that increasing η/s to 2/4π can lower v_2 and v_3 for all p_T, though the effect is stronger in higher p_T, which also agrees with hydrodynamic approach.



It is interesting to study the difference between the effect of η/s on collective flows for small colliding systems and large ones.
To this end, we show in Fig. <ref> the ratio v_n{2} (4πη/s=2)/v_n{2} (4πη/s=1) as a function of transverse momentum in 5.02 TeV p+Pb and 5.02 TeV Pb+Pb collisions at 0-1% centrality with same parameters. We find that in general increasing η/s leads to the relative smaller effect on the decrease of v_2 and relative larger effect on the decrease of v_3 for all p_T. However, it is found that this decrease is stronger at larger p_T in small colliding systems, while it is almost uniform for all p_T in large colliding systems. We find that the effect of η/s is stronger for p+Pb collisions, which can be seen by the solid red and blue lines in Fig. <ref>.
In fact the effect in Pb+Pb is about a 10% for v_2 (dashed red line) and 20% for v_3 (dahsed blue line)
while for p+Pb is about twice as large already at p_T ≃ 1 GeV/c and further increase with p_T.
This suggest that ultra-central p+Pb collisions may supply more sensitivity for the determination of the η/s
w.r.t. to Pb+Pb collisions. However this larger sensitivity has been seen in our calculation only for the most central 
collisions, while at large p+Pb centrality above 20-30 % (small multiplicity N_ch< 80) the impact
of η/s becomes similar to Pb+Pb collisions if not smaller.



To this end we have to add that the initial eccentricities in p+Pb collisions are strongly dependent on the width σ
of the fluctuation of transverse density, as it is shown in the upper panel of Fig. <ref>.
It is seen that there are dramatic enhancements of ⟨ϵ_2⟩ and ⟨ϵ_3⟩ with smaller σ.
Because v_n and ϵ_n are strong correlated, we notice that in principle an increasing η/s can be compensated
by a smaller σ width of the spacial fluctuations that induce a larger ϵ_n. However it is interesting to mention that a flat behavior of ⟨ϵ_2⟩≈0.28 with centrality is in nice agreement with the behavior extracted from a recent analysis by kinetic theory of pA collisions presented in Ref. <cit.>. 

In the lower panel of Fig. <ref>, we also show ϵ_2 and ϵ_3 as a function of centrality in Pb+Pb collisions with σ=0.55 and 0.25 fm; unlikely to p+Pb collisions, there is almost no dependence of eccentricity on σ. Therefore while in AA collisions one has a very limited dependence on the width of the initial state fluctuations this
is large in pA collisions. This limits the possibility to costrain the η/s of the matter created if there are no independent
observables that allow to constraint also the initial state fluctuations. In Fig. <ref> we show the ratio of the anisotropies  as a function of transverse momentum
v_n{2} (σ_1)/v_n{2} (σ_2) when changing the widths from a σ_1=0.55 fm
to σ_2=0.45 fm in ultra-central collisions(0-0.3%), assuming an η/s=1/4π. 
The main information we get is that the correction we have is essentially p_T independent especially for v_2
and this shows that changing σ is not equivalent to change the η/s that induce a quite p_T
dependent modification of v_2(p_T). Therefore the two effect does not compensate each other when
looking at the p_T dependence of the anisotropies.



We also report in Fig. <ref> that doing a study of v_2 and v_3 as a function of the multiplicity
 of charged hadrons at |η|<2.4 and 0.3<p_T<3 GeV/c
 we find a good agreement also for the integrated v_n, but  in the most central collisions while we tend to underestimate v_2 in more peripheral collisions.
We find that both v_2 and v_3 increase with the increase of charged particle multiplicity. Because the eccentricity and size of the p+Pb systems change less than 6% from 0-0.3% to 30-40 % centrality, the only reason for this is the shorter lifetime of low multiplicity events due to a lower initial energy density. 


However a small change of  σ from 0.55 fm to 0.45 fm, allows to reproduce the v_2 in low multiplicity events, while we overestimate it in high multiplicity events. Changing σ from 0.55 fm to 0.45 fm doesn't change 
significantly the integrated v_3, which can be seen by the relative small differences between solid and dashed blue lines in Fig. <ref>.
A change in σ as a function of centrality does not have a solid physical motivation, we only wanted to convey the message that the quantitative agreement can be to some extent
tuned by the size of the spatial fluctuations. 
A real understanding of pA asks in the future at least to extend the study to higher harmonics and their correlations as done
for AA collisions.









§ THE IMPACT OF MINIJET AT MID-P_T


In the above section, we find that our results on spectra of charged hadrons at higher p_T 
significantly underestimate experimental measurements, see Fig. <ref>. 
However a main advantage of the transport approach with respect to the pure hydrodynamics is the possibility
to self-consistently include initial conditions that significantly deviate from the equilibrium one as the power law tail
associated to mini-jet production. In this section we extend our study including minijet at moderate mid transverse momentum,
however still limiting ourself to a region of p_T < 5 GeV/c because of limitations due to statistics but also because at higher
p_T would be necessary to properly include also the radiative energy loss and the details that distinguish it
from an elastic collisional energy loss.
Because the non-equilibrium production of minijets in initial conditions can have an impact on both spectra and collective flows at higher p_T, we thus try to include them in initial conditions in this section. To do this, we use the spectra of gluons and quarks from CUJET Collaboration <cit.>  in 5.02 TeV p+p collisions, which is parametrized as 

    dN_g,q/d^2p_Tdy=1/120(a/1+p_T/b)^c (GeV/c)^-2,

where a=3.36479, b=1.53518  GeV/c, c=6.16767 for gluons and a=2.24291, b=1.31444 GeV/c, c=5.72108 for three flavors light quarks (same for anti-quarks). The production of minijets in 5.02 TeV p+Pb collisions should be scaled with binary collisions, which can be obtained from the Monte Carlo Glauber model described above. For the distribution of minijets in transverse plane, we assume they are centered around the center of each binary collision pair, with the same Gaussian width as thermal partons. On the other hand, we assume minijets have no asymmetric distributions in spatial rapidity as they are produced by binary collisions. For the longitudinal momentum distribution, we assume that their rapidity is same as spatial rapidity because they are not in thermal equilibrium.

We have found that including minijets in initial conditions can lead to the enhancement of spectra at higher p_T, which can be seen by the solid black line in Fig. <ref>, where we include minijets at p_T>3 GeV/c. The reason for this is trivial, because exponential law decay of thermal partons decays faster compared to the power law decay of minijets.


Moreover the mini-jets not only enhance the spectra of charged hadrons at higher p_T, but also can affect the transverse momentum dependence of collective flows. It is seen in Fig. <ref> by the solid red and blue lines that the inclusion of minijets decreases v_2 and v_3 at p_T ≃ 2.5 GeV/c, which leads to a better agreement with experimental data up to 5 GeV/c.

 
To our knowledge this is the first time the collective flows v_2,3 in pA collisions are predicted correctly in such a
wide range of p_T. We think that the agreement with the experimental data shown in Fig. <ref> 
further strengthen the validity of the interpretation of the anisotropies observed as coming from an hydro-like expansion. We have to clarify that by hydro-like behavior we simply mean that we can achieve a reasonable description of collective anisotropies only assuming a matter expanding with a significantly highly rate of collisions similarly to a fluid. However, it remains an important open question if such an expanding phase is fully hydrodynamical or it is consistent with a transport evolution where non-hydrodynamical (particle-like) collective excitations are relevant. Very recently in Ref. <cit.> a new quantitatively analysis, employing a conformal kinetic theory, has shown that pA may lie at least in a region of transition between pure hydrodynamics and kinetic-particle evolution. It will certain be a study to be pursued to extend the study to our approach that contains also non conformal contribution to the expansion and could investigate also potential impact of the power-law tail in the distribution function.




§ FLOW CORRELATIONS BETWEEN V_N AND Ε_N FROM LARGE COLLIDING SYSTEMS TO SMALL ONES

The correlations between collective flows v_n and initial asymmetry in coordinate space ϵ_n 
in AA collisions have been studied in event-by-event hydrodynamics and transport approaches in recent years <cit.>. In general it has been shown that v_2 is strongly correlated with ϵ_2 while higher flows are correlated weaker with ϵ_n>2. In this section, we will 
study, to our knowledge for the first time in pA collisions,  
compare such correlations to the one in AA collisions.

To characterize the strength of the correlation between v_n and ϵ_n, we adopt the same correlation coefficient C(n,m) defined in Ref. <cit.>

    C(n,m)=∑_i(ϵ_n^i-⟨ϵ_n⟩)(v_m^i-⟨ v_m⟩)/√(∑_i(ϵ_n^i-⟨ϵ_n⟩)^2∑_i(v_m^i-⟨ v_m⟩)^2),

where ϵ_n is defined as

    ϵ_n=√(⟨ r_T^ncos (nϕ)⟩^2+⟨ r_T^n sin(nϕ)⟩^2)/⟨ r_T^n ⟩.

C(n,m) close to one corresponds to the stronger linear correlation between initial ϵ_n and final v_m.


In Fig. <ref> we show the correlations between v_2, v_3 and ϵ_2 and ϵ_3 for 5.02 TeV p+Pb collisions at 0-0.3% and 20-30% centrality class with Gaussian width σ=0.55 fm. In the upper panel of Fig. <ref>, it is seen that the correlation between ϵ_2 and v_2 is larger than the one between ϵ_3 and v_3. From the central collision to peripheral collision, we find that the correlation decreases both for n=2 and n=3, though it decreases faster for the correlation between ϵ_3 and v_3. This is a little different from the trend in Pb+Pb collisions, where there is almost no drop of the correlation for n=2, and only a small drop of it for n=3 <cit.>. In Fig. <ref>, we furthermore indicate the values of ⟨ v_n⟩/⟨ϵ_n⟩ ratio, which decreases also for more peripheral collisions, same as the correlations.


We also studied the effect of varying the gaussian width σ to see the effect on of  C(n,n) and ⟨ v_n⟩/⟨ϵ_n⟩ ratio in p+Pb collisions, and it is found that the their values are very similar for 
both the centrality class  0-0.3% and 30-40%. Because of this, we only compare our results in p+Pb collisions with σ=0.55 fm to what calculated in Ref. <cit.> in Pb+Pb collisions.



In Fig. <ref>, we show the ⟨ v_n⟩/⟨ϵ_n⟩ ratio as a function of centrality class in 5.02 TeV p+Pb collisions for n=2 and 3 as well as those calculated  in Ref. <cit.> in 5.02 TeV for Pb+Pb collisions. It is seen that both  ⟨ v_2⟩/⟨ϵ_2⟩ and ⟨ v_3⟩/⟨ϵ_3⟩ are quite smaller in p+Pb collisions compared to Pb+Pb collisions. Furthermore they all decrease at increasing of centrality, regardless of the size of colliding systems, but they decreases quite faster in small colliding systems.
From Fig. <ref> we can see that the smaller v_2,3 in p+Pb is due to a quite reduced efficiency in converting
the initial eccentricities ϵ_n even if the η/s assumed is the same as the one for Pb+Pb. We notice that the value of ⟨ v_2⟩/⟨ϵ_2⟩ and ⟨ v_3⟩/⟨ϵ_3⟩, representing the efficiency of conversions of space anisotropy, are quite similar to those obtained in viscous hydrodynamics in Ref. <cit.> that however is a study at RHIC energy for p+Au, d+Au and ^3He+Au.


In order to better visualize the relation between the correlation coefficient C(n,n) and centrality, we plot in Fig. <ref> the correlation coefficient as a function of centrality in 5.02 TeV p+Pb collisions for n=2, 3. It is seen by the solid red and blue lines that both C(2,2) and C(3,3) decreases with increase of centrality, while the correlation for triangular flow C(3,3) is smaller than elliptic flow C(2,2), and decreases faster. Compared to Pb+Pb collisions from Ref. <cit.>, we find that the correlation coefficients for n=2 and n=3 are smaller for small colliding systems and decreases faster with centrality class.


We propose to access such correlations in pA systems because their measurements could give an important contribution to
validate the interpretation of the anisotropies v_n in pA collisions as hydro-like collective expansion.

The correlation between collective flows and the initial geometry is strong, and the initial geometry of p+Pb collisions is dominated by fluctuations, the distribution of v_n should be quite relevant as an indicator of initial conditions. 
Therefore we conclude our study of v_2,3 in pA collisions discussing the normalized variance
 σ_v_n/⟨ v_n⟩. In Fig. <ref> 
we plot the centrality dependence of  both σ_v_n/⟨ v_n⟩ and σ_ϵ_n/⟨ϵ_n⟩ in 5.02 TeV p+Pb collisions, as well as the same one in 5.02 TeV Pb+Pb collisions from Ref. <cit.>. Shown in solid and dashed blue lines in the upper panel of Fig. <ref> for n=2, we find that σ_ϵ_2/⟨ϵ_2⟩ increases slightly with centrality class in p+Pb collisions. 
In Pb+Pb the decrease with centrality class can be attributed to 
the additional contribution of ϵ_2 from the global average geometry in latter case. 
According to our modeling on intial conditions in pA ϵ_n is dominated by fluctuations at all centralities
entailing a σ_v_2/⟨ v_2⟩ that even slightly increase with centrality 
as shown by the solid red line, compared to the decrease of it with centrality in large colliding systems 
as shown by the dashed red line.
Of course observing such a different trend would give further support to the current modeling of the initial conditions.



In the lower panel of Fig. <ref>, we show σ_v_3/⟨ v_3⟩ and σ_ϵ_3/⟨ϵ_3⟩ as a function of centrality class in p+Pb and Pb+Pb collisions. It is seen by the solid and dashed red lines that σ_v_3/⟨ v_3⟩ increases with centrality class in p+Pb collisions while stays almost constant in Pb+Pb collisions, which indicates that triangular flow is dominated by fluctuations in both small and large colliding systems, but
anyway there is significant tendency to stay larger in pA systems. We also find that  σ_ϵ_3/⟨ϵ_3⟩ is smaller in large colliding systems compared to small ones. Shown by the solid and dashed red lines, σ_v_3/⟨ v_3⟩ also increases slowly with centrality class in p+Pb collisions, while it stays the same in Pb+Pb collisions and below the value in p+Pb collisions, which agrees with the trend of σ_ϵ_3/⟨ϵ_3⟩.
We think that such a measurements can significantly contribute to validate or falsify our modeling of
the initial conditions and of the dynamics of pA collisions.



§ CONCLUSIONS AND DISCUSSIONS

We use an event-by-event transport approach whose initial conditions are generated by wounded quark model as well as negative binomial distribution overlaid on partons production for each participant quark, to study several hadron observables in p+Pb collisions. We find that we can reproduce quite well dN/dη, the distribution of charged particle multiplicity and the spectra of charged hadrons up to 1.5 GeV/c at midrapidity. 
We have shown that the anisotropies v_2 and v_3 are building-up in time in a way that roughly scales with
the mean square radius of the system. So the formation time in pA system will be in general about a factor 4 faster
with respect to most central AA collisions.
A main result is that the transverse momentum dependence of v_2 and v_3 predicted in our approach
agree with experimental measurements with η/s=1/4π up to about 2.5 GeV/c, as well as the integrated v_2 and v_3. 
This is in general confirm several results already obtained in viscous hydrodynamics <cit.>, we also find that however ultra-central p+Pb collisions are more sensitive to
the value of the η/s assumed.
A specific advantage of our approach is the possibility to include also power law tail associated to mini-jets production,
we have found that this allows to extend the 
 agreement with experimental data on  collective flows at higher p_T at least up to
 about 5 GeV/c and this is also accompanied by a much better prediction of the transverse momentum spectrum.
We have also pointed out that at variance with respect to Pb+Pb in small systems the ϵ_n and v_n
are strongly dependent on the width of the initial spatial gaussian fluctuations. There is some interplay between
the value of the η/s and the width of the initial fluctuations, however especially for v_2 reducing the width
induces an increase of the elliptic flow that is p_T independent.

Finally, we furthermore study the correlations between v_n and ϵ_n that has been intensively studied for
AA collisions, but to our knowledge they are for the first time discussed for pA collisions. 
In general, we find that the correlation coefficient is still strong for n=2 in 0-0.3% collisions similarly
to AA collisions, but it becomes weaker for either the higher harmonics v_3 or larger centrality class, and the correlation coefficient decreases quite faster with centrality in p+Pb collisions compared to Pb+Pb collisions. So we can say that except
ultra-central collisions and only for v_2 we expect significantly less correlations in p+Pb collisions.
Besides that, we further predict that except for ultra-central collisions the variance σ_v_2/⟨ v_2⟩ should be nearly centrality independent and quite larger than the one in Pb+Pb.
We plan in a further study to investigate if an initial glasma phase can affect such behaviors. In the meantime
experimental measurements can shed new light on the goodness of the present modeling of pA as liquid drops
with local density fluctuations expanding nearly hydrodynamically.



§ ACKNOWLEDGEMENTS

We gratefully acknowledge useful discussions with Jamie L. Nagle. The work of  Y.S. is supported by a INFN post-doc fellowship within the national SIM project.




