Study of collective anisotropies $v_2$ and $v_3$ and their fluctuations in $pA$ collisions at LHC within a relativistic transport approach

TL;DR

This study develops a relativistic transport model incorporating initial fluctuations to analyze collective flow anisotropies in p+Pb collisions at LHC, successfully matching experimental data and exploring the effects of shear viscosity and initial conditions.

Contribution

The paper introduces a relativistic transport approach with initial fluctuations and minijet contributions, extending the understanding of flow anisotropies in small collision systems.

Findings

The model accurately predicts experimental flow measurements.

Including minijet effects improves high p_T data agreement.

Flow fluctuation patterns differ between p+Pb and Pb+Pb collisions.

Abstract

We have developed a relativistic transport approach at fixed $\eta/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 $\eta/s=1/4\pi$, a result similar to previous studies within a viscous hydrodynamics approach. Besides, we further discuss the sensitivity of the results on both $\eta/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 $\sigma_{v_n}/\langle v_n\rangle$ has a quite different evolution with centrality for p+Pb, so their measurement could provide some further hint about the correctness of current modelling.

Equations

$$p_{\mu }{\partial }^{\mu }f(x,p)+m^{\ast }{\partial }^{\mu }{m}^{\ast }{\partial }_{p^{\mu }}f(x,p)=\mathcal{C}[f]$$

$$\eta =f(z)\frac{T}{{\sigma }_{tot}}$$

$$f(z)=\frac{15}{16}\frac{{\left[z^2{K}_3(z)\right]}^2}{(15z^2+2)K_2(2z)+(3z^3+49z)K_3(2z)}$$

$${\rho }_{L\pm }(\eta )=(1\pm \frac{\eta }{{\eta }_m})\exp (-\frac{(|\eta |-{\eta }_0{)}^2}{2{\sigma }_{\eta }^2}\theta (|\eta |-{\eta }_0)),$$

$$\frac{dN}{d^2{\mathbf{x}}_{\perp }d\eta }=\sum _{i=1}^{N_{\mathrm{part}}}{n}_i{\rho }_{\perp }({\mathbf{x}}_{\perp }-{\mathbf{x}}_{\mathbf{i}}){\rho }_{L\pm }(\eta ),$$

$$\frac{dN}{d^2{\mathbf{x}}_{\perp }d\eta }=\frac{g{\tau }_0}{(2\pi {)}^2}{p}_T{m}_T$$

$$\times \cosh (\eta -\mathrm{y}){\mathrm{e}}^{-{\mathrm{m}}_{\mathrm{T}}\cosh (\eta -\mathrm{y})/T({\mathbf{x}}_{\perp },\eta )}\mathrm{d}{\mathrm{p}}_{\mathrm{T}}\mathrm{dy},$$

$$\frac{d{N}_{g,q}}{d^2{p}_{T}dy}=\frac{1}{120}(\frac{a}{1+p_T/b}{)}^c(\mathrm{GeV}/c{)}^{-2},$$

$$C(n,m)=\frac{{\sum }_i({ϵ}_n^i-⟨{ϵ}_n⟩)(v_m^i-⟨v_m⟩)}{\sqrt{{\sum }_i({ϵ}_n^i-⟨{ϵ}_n⟩{)}^2{\sum }_i(v_m^i-⟨v_m⟩{)}^2}},$$

$${ϵ}_n=\frac{\sqrt{⟨r_T^n\cos (\mathrm{n}\varphi ){⟩}^2+⟨{\mathrm{r}}_{\mathrm{T}}^{\mathrm{n}}\sin (\mathrm{n}\varphi ){⟩}^2}}{⟨r_T^n⟩}.$$

Full text

t

collectivity, flow fluctuations, small colliding systems

I Introduction

High energy nucleus-nucleus collision experiments at the BNL Relativistic Heavy Ion Collider (RHIC)  Adams et al. (2005); Adcox et al. (2005); Arsene et al. (2005); Back et al. (2005) and the CERN Large Hadron Collider (LHC) Aamodt et al. (2008) 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 Romatschke and Romatschke (2007); Song and Heinz (2008); Schenke et al. (2011); Gale et al. (2013); Song et al. (2011); Niemi et al. (2016) and transport approach Ferini et al. (2009); Xu and Greiner (2009); Ruggieri et al. (2014); Plumari et al. (2015a); Konchakovski et al. (2012); Karpenko et al. (2013) have shown that such behavior is consistent with a matter with low shear viscosity to specific entropy ratio $\eta/s\simeq\,0.1-0.2$ close to the conjectured lower bound $\eta/s=1/4\pi$ for a strongly interacting system Kovtun et al. (2005).

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 Chatrchyan et al. (2013a); Abelev et al. (2013); Aad et al. (2013); Adare et al. (2013). Several theoretical studies based on hydrodynamics Bozek (2012); Bozek and Broniowski (2013); Kozlov et al. (2014); Shen et al. (2017); Weller and Romatschke (2017) and transport approach Bzdak and Ma (2014); Greif et al. (2017) 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 $\eta/s$. Such an approach has been shown to reproduce the viscous hydrodynamics system expansion in $AA$ collisions Ruggieri et al. (2014); Plumari et al. (2015a) and even the ideal hydrodynamics behavior in the infinite cross section limit Plumari et al. (2015b). 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 $\eta/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 McLerran and Skokov (2017); Mace et al. (2018); Greif et al. (2017). 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.

II Model setup

To model the evolution of the small collision systems, we are employing the same relativistic transport code developed in these years Ferini et al. (2009); Plumari and Greco (2012); Plumari et al. (2012); Ruggieri et al. (2013, 2014); Plumari et al. (2015a), 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 Xu and Greiner (2005); Ferini et al. (2009); Plumari et al. (2012).

[TABLE]

The key aspect of our approach is to gauge the collision integral ${\cal C}[f(x,p)]$ locally to a specific value of the viscosity to entropy density ratio $\eta/s$  Plumari et al. (2012). Furthermore we can switch smoothly it to an increasing $\eta/s$ one, reaching the estimated value in hadronic phase Chen et al. (2007); Demir and Bass (2009), when the system reaches the cross over region. This is realized by determining the total isotropic cross section according to the Chapmann-Enskog approximation:

[TABLE]

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

[TABLE]

where $K_{n}$-s are the modified Bessel functions. The entropy density for a massive system is given by s=\rho\big{(}4+zK_{1}(z)/K_{2}(z)\big{)}. In the ultra-relativistic limit $z\to 0$ and the function $f(z)\to 1.2$ and we recover the massless limit for the $\eta$. As shown in Plumari et al. (2012) The expression for $\eta$ in Eq.(2) is in quite good agreement at level of $3-5\%$ with the Green-Kubo formula Plumari et al. (2012). 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 Bialas and Czyz (1979); Anisovich et al. (1978); Eremin and Voloshin (2003); Bozek et al. (2016), 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=\frac{r^{2}}{r_{0}^{3}}e^{-r/r_{0}}$ with $r_{0}=0.3$ fm Bozek et al. (2016), 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^{-\pi r^{2}/\sigma_{qq}}$ with $\sigma_{qq}=13.6$ mb in 5.02 TeV p+Pb collision Bozek et al. (2016).

For the distribution of the spatial rapidity, we take the profile from Ref. Bozek and Wyskiel (2010)

[TABLE]

where the parameters are chosen as $\eta_{m}=5.7$, $\eta_{0}=2.5$ and $\sigma_{\eta}=2.5$ to have the same shape of $dN_{\rm{ch}}/d\eta$ measured by ATLAS Collaboration Milov (2014) as shown in Fig. 1. 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

[TABLE]

where $\mathbf{x_{i}}$ is the transverse position of each participant quark, $n_{i}$ is the number of partons generated by each participant, and $\rho_{\perp}(\mathbf{x_{\perp}})=\frac{1}{2\pi\sigma^{2}}e^{-\frac{\mathbf{x_{\perp}}^{2}}{2\sigma^{2}}}$. For the Gaussian distribution of the parton transverse density, we will change the parameter $\sigma$ 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. (5) to be $n_{0}N$, where $N$ is sampled according to NBD $P(N)=\frac{\Gamma(N+\kappa){\overline{n}}^{N}\kappa^{\kappa}}{\Gamma(\kappa)N!(\overline{n}+\kappa)^{N+\kappa}}$, and $n_{0}$ is a constant such that the final charged particle multiplicity is same as that measured in experiments. We find that $\kappa=0.54$, $\overline{n}=3.9$ and $n_{0}\approx 2.352$ can almost reproduce the distribution of charged particle measured by CMS Collaboration Chatrchyan et al. (2013b), as shown in Fig. 2. We also mention that our $N_{\rm{ch}}$ is assumed to be equal to $N_{\rm{trk}}^{\rm{offline}}$ after the efficiency corrections usually adopted in hydro calculation Kozlov et al. (2014); Shen et al. (2017); Bozek et al. (2016). 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 Bozek et al. (2016).

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

[TABLE]

where $g=2\times 8+3\times 2\times 6=52$ is the degree of freedom of partons (three flavor quarks and gluons), $\tau_{0}=0.4$ fm$/c$ is the thermalization time (again taken as in standard hydro approach for $pA$), and $m_{T}=\sqrt{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. 1. Also in thermal equilibrium it generates an equation of state close to the one calculated in lattice QCD Plumari (2019). We use the above relation to calculate the temperature locally from initial parton density, and then sample the momenta of partons according to Eq. (6).

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.

III 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_{\rm test}$ appropriate to compute observables, where $N_{\rm 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_{\rm{test}}=4000$ and $\Delta x=\Delta y=0.15$ fm in the simulations guarantee a convergency for $v_{2,3}(p_{T})$ up to $p_{T}\sim$5 GeV$/c$.

In Fig. 3, we compare our spectra to experimental measurements in minimum-bias 5.02 TeV p+Pb collisions from ALICE Collaboration Abelev et al. (2014), 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 Bozek (2012); Shen et al. (2017). 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.10.

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 Bhalerao et al. (2005) 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 $\rm R=\sqrt{\langle r^{2}\rangle}$, where $\langle r^{2}\rangle$ is the mean square radius of the system weighted with the local density of the system. Indeed the value of $\rm 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. 4 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 $\sigma=0.55$ fm for gaussian fluctuations of the parton transverse density and $\eta/s=1/4\pi$, we find also a good agreement with experimental data measured by CMS Collaboration Chatrchyan et al. (2013b) for elliptic and triangular flows, as shown in Fig. 5. It is seen by the solid red and blue lines that $v_{2}$ and $v_{3}$ agrees with experimental data up to $p_{T}\simeq 2.5$ GeV$/c$ in our calculations which is quite similar to the results obtained within viscous hydrodynamics with very similar initial conditions Bozek (2012); Bozek and Broniowski (2013); Shen et al. (2017). 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 $\eta/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 $\eta/s$ on $v_{2}$ and $v_{3}$, which is shown by the dashed red and blue lines in Fig. 5. It is seen that increasing $\eta/s$ to $2/4\pi$ 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 $\eta/s$ on collective flows for small colliding systems and large ones. To this end, we show in Fig. 6 the ratio $v_{n}\{2\}(4\pi\eta/s=2)/v_{n}\{2\}(4\pi\eta/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 $\eta/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 $\eta/s$ is stronger for p+Pb collisions, which can be seen by the solid red and blue lines in Fig. 6. 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}\simeq$ 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 $\eta/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_{\rm{ch}}<80$) the impact of $\eta/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 $\sigma$ of the fluctuation of transverse density, as it is shown in the upper panel of Fig. 7. It is seen that there are dramatic enhancements of $\langle\epsilon_{2}\rangle$ and $\langle\epsilon_{3}\rangle$ with smaller $\sigma$. Because $v_{n}$ and $\epsilon_{n}$ are strong correlated, we notice that in principle an increasing $\eta/s$ can be compensated by a smaller $\sigma$ width of the spacial fluctuations that induce a larger $\epsilon_{n}$. However it is interesting to mention that a flat behavior of $\langle\epsilon_{2}\rangle\approx 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. Kurkela et al. (2019).

In the lower panel of Fig. 7, we also show $\epsilon_{2}$ and $\epsilon_{3}$ as a function of centrality in Pb+Pb collisions with $\sigma=0.55$ and $0.25$ fm; unlikely to p+Pb collisions, there is almost no dependence of eccentricity on $\sigma$. 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 $\eta/s$ of the matter created if there are no independent observables that allow to constraint also the initial state fluctuations. In Fig. 8 we show the ratio of the anisotropies as a function of transverse momentum $v_{n}\{2\}(\sigma_{1})/v_{n}\{2\}(\sigma_{2})$ when changing the widths from a $\sigma_{1}=0.55$ fm to $\sigma_{2}=0.45$ fm in ultra-central collisions(0-0.3%), assuming an $\eta/s=1/4\pi$. 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 $\sigma$ is not equivalent to change the $\eta/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. 9 that doing a study of $v_{2}$ and $v_{3}$ as a function of the multiplicity of charged hadrons at $|\eta|<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 $\sigma$ 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 $\sigma$ 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. 9. A change in $\sigma$ 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.

IV 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. 3. 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 Xu et al. (2016) in 5.02 TeV p+p collisions, which is parametrized as

[TABLE]

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. 10, 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. 11 by the solid red and blue lines that the inclusion of minijets decreases $v_{2}$ and $v_{3}$ at $p_{T}\simeq$ 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. 11 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. Kurkela et al. (2019) 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.

V Flow correlations between $v_{n}$ and $\epsilon_{n}$ from large colliding systems to small ones

The correlations between collective flows $v_{n}$ and initial asymmetry in coordinate space $\epsilon_{n}$ in $AA$ collisions have been studied in event-by-event hydrodynamics and transport approaches in recent years Gardim et al. (2012); Chaudhuri et al. (2013); Niemi et al. (2013); Plumari et al. (2015a). In general it has been shown that $v_{2}$ is strongly correlated with $\epsilon_{2}$ while higher flows are correlated weaker with $\epsilon_{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 $\epsilon_{n}$, we adopt the same correlation coefficient $C(n,m)$ defined in Ref. Plumari et al. (2015a)

[TABLE]

where $\epsilon_{n}$ is defined as

[TABLE]

$C(n,m)$ close to one corresponds to the stronger linear correlation between initial $\epsilon_{n}$ and final $v_{m}$.

In Fig. 12 we show the correlations between $v_{2}$, $v_{3}$ and $\epsilon_{2}$ and $\epsilon_{3}$ for 5.02 TeV p+Pb collisions at 0-0.3% and 20-30% centrality class with Gaussian width $\sigma=0.55$ fm. In the upper panel of Fig. 12, it is seen that the correlation between $\epsilon_{2}$ and $v_{2}$ is larger than the one between $\epsilon_{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 $\epsilon_{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$ Plumari et al. (2015a). In Fig. 12, we furthermore indicate the values of $\langle v_{n}\rangle/\langle\epsilon_{n}\rangle$ ratio, which decreases also for more peripheral collisions, same as the correlations.

We also studied the effect of varying the gaussian width $\sigma$ to see the effect on of $C(n,n)$ and $\langle v_{n}\rangle/\langle\epsilon_{n}\rangle$ 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 $\sigma=0.55$ fm to what calculated in Ref. Plumari et al. (2015a) in Pb+Pb collisions.

In Fig. 13, we show the $\langle v_{n}\rangle/\langle\epsilon_{n}\rangle$ 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. Plumari et al. (2015a) in 5.02 TeV for Pb+Pb collisions. It is seen that both $\langle v_{2}\rangle/\langle\epsilon_{2}\rangle$ and $\langle v_{3}\rangle/\langle\epsilon_{3}\rangle$ 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. 13 we can see that the smaller $v_{2,3}$ in p+Pb is due to a quite reduced efficiency in converting the initial eccentricities $\epsilon_{n}$ even if the $\eta/s$ assumed is the same as the one for Pb+Pb. We notice that the value of $\langle v_{2}\rangle/\langle\epsilon_{2}\rangle$ and $\langle v_{3}\rangle/\langle\epsilon_{3}\rangle$, representing the efficiency of conversions of space anisotropy, are quite similar to those obtained in viscous hydrodynamics in Ref. Nagle et al. (2014) 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. 14 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. Plumari et al. (2015a), 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 $\sigma_{v_{n}}/\langle v_{n}\rangle$. In Fig. 15 we plot the centrality dependence of both $\sigma_{v_{n}}/\langle v_{n}\rangle$ and $\sigma_{\epsilon_{n}}/\langle\epsilon_{n}\rangle$ in 5.02 TeV p+Pb collisions, as well as the same one in 5.02 TeV Pb+Pb collisions from Ref. Plumari et al. (2015a). Shown in solid and dashed blue lines in the upper panel of Fig. 15 for $n=2$, we find that $\sigma_{\epsilon_{2}}/\langle\epsilon_{2}\rangle$ 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 $\epsilon_{2}$ from the global average geometry in latter case. According to our modeling on intial conditions in $pA$ $\epsilon_{n}$ is dominated by fluctuations at all centralities entailing a $\sigma_{v_{2}}/\langle v_{2}\rangle$ 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. 15, we show $\sigma_{v_{3}}/\langle v_{3}\rangle$ and $\sigma_{\epsilon_{3}}/\langle\epsilon_{3}\rangle$ as a function of centrality class in p+Pb and Pb+Pb collisions. It is seen by the solid and dashed red lines that $\sigma_{v_{3}}/\langle v_{3}\rangle$ 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 $\sigma_{\epsilon_{3}}/\langle\epsilon_{3}\rangle$ is smaller in large colliding systems compared to small ones. Shown by the solid and dashed red lines, $\sigma_{v_{3}}/\langle v_{3}\rangle$ 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 $\sigma_{\epsilon_{3}}/\langle\epsilon_{3}\rangle$. 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.

VI 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\eta$, 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 $\eta/s=1/4\pi$ 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 Bozek (2012); Bozek and Broniowski (2013); Kozlov et al. (2014); Weller and Romatschke (2017); Shen et al. (2017), we also find that however ultra-central p+Pb collisions are more sensitive to the value of the $\eta/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 $\epsilon_{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 $\eta/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 $\epsilon_{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 $\sigma_{v_{2}}/\langle v_{2}\rangle$ 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.