Glauber model for small system using anisotropic and inhomogeneous density profile of proton

Recent studies reveal that at high energies, collisions of small system like p + p give signatures similar to that widely observed in heavy ion collisions hinting towards a possibility of forming a medium with collective behaviour. With this motivation, in this work, we have used Glauber model, which is traditionally applied to heavy ion collisions, in small system using anisotropic and inhomogeneous density profile of proton and found that the proposed model reproduces the charged particle multiplicity distribution of p + p collisions at LHC energies very well. Collision geometric properties like mean impact parameter, mean number of binary collisions (〈Ncoll〉) and mean number of participants (〈Npart〉) at different multiplicities are determined. Having estimated 〈Ncoll〉, we have calculated nuclear modification-like factor (Rpp) in p+p collisions. We also estimated eccentricity and elliptic flow as a function of charged particle multiplicity using linear response to initial geometry.


I. INTRODUCTION
Results of relativistic proton proton (p + p) collisions are used as reference or base line for interpreting various results of heavy ion collisions at relativistic energies, which are aimed at creation and characterization of phases of strongly interacting matter governed by Quantum Chromodynamics (QCD). At high temperature/density, a deconfined thermalised state of quarks and gluons called Quarkgluon plasma (QGP) has been predicated by lattice QCD based calculations [1,2]. Values of ratio of certain observable measured in heavy ion collisions, such as number of produced strange particles, production of J/ψ, to those measured in p + p collision are interpreted as signature of partonic medium formation in heavy ion collisions; e.g., enhancement of number of strange particles and suppression in number of J/ψ in collisions of heavy ion with respect to that of p + p are taken as signatures of QGP formation in relativistic heavy collisions [3][4][5][6][7][8][9][10][11]. Apart from taking values of such ratio as confirmation for creation of QGP in such collisions, they are also used in characterizing QGP as well as in verifying and constraining different theoretical models. For such interpretations, it is assumed that, in p + p collisions, no intermediate partonic medium is formed. However, recent results show that such assumptions * Electronic address: Raghunath.Sahoo@cern.ch † Electronic address: jane@vecc.gov.in may not be correct for high multiplicity p + p collisions [12][13][14]. In this regard, for understanding dynamics of QCD medium formed in relativistic heavy ion collisions, it is very crucial to understand results of p + p collisions.
For understanding underlying physics, and for characterizing the possibly formed medium in these collisions by using theoretical models, proper estimation of initial condition for such collisions is crucially important. In the context of heavy ion collisions, the use of initial condition is very crucial for extracting all the signals of QGP [15][16][17]. For example, the elliptic flow (v 2 ) of hadrons are calculated using the transverse momentum (p T ) spectrum of hadrons. The p T spectra are estimated using the hydrodynamical models. To solve the hydrodynamical equations, initial conditions and equation of state are required as inputs. Therefore, the shear viscosity of the system extracted through v 2 will be sensitive to initial conditions, i.e., any uncertainties in initial condition will be reflected in shear viscosity [18].
In high energy heavy-ion collisions, interpretation of results relies on the use of a model based on initial matter distribution resulting from the overlap of the two colliding nuclei at a given impact parameter (b). Indeed, for estimating quantities such as: (i) the centrality dependence of any observable expressed by the number of participating nucleons in the collision, N part (b), (ii) the number of binary nucleon-nucleon collisions, N coll (b) used to derive the nuclear modification factor (R AA ) from the ratio of AA over pp spectra, (iii) the elliptic and triangular flow parameters (v 2 ) and (v 3 ) normalized by the eccentricity 2 (b) and triangulation 3 (b) of the overlap region, and (iv) the average surface area, A(b) and (v) path length, L(b) of the interaction region, knowing the nuclear overlap function T AA (b) is important. And, this overlap function depends on a realistic model of the collision geometry [19].
In a similar way, to understand, one of the most recent on going discussed topic, the possible medium formation in high-multiplicity p+p collisions (events) in light of fluid dynamic models, knowing proper initial conditions is important. Apart from this, knowing proper initial condition can also give a possible way to define centrality classes and the base needed for properly defining suppression factors or ratios for comparing results of event of different multiplicity classes produced in p+p collisions [20]. Appropriate initial conditions can be chosen by considering that it should reproduce certain aspects of results such as multiplicity distribution or centrality distribution of various observables related to the events.
For constructing proper initial conditions for p+p collisions, at first attempt, one follows the way similar to that of heavy ion collisions. Initial conditions for heavy ion collisions are modeled in two kinds of distinct approaches-in one, one considers nucleonic or partonic collisions for energy deposition in the collision zone, those are based on Glauber model [31], and in second kinds of approach QCD based calculations are employed to estimate initial energy deposition by gluonic fields originated from partonic currents of colliding nuclei. So these will also be obvious approaches for modeling initial conditions in p+p collisions. As models based on Glauber modeling, are very successful in reproducing various results of relativistic heavy ion collisions, one can consider models for initial conditions of p + p collisions which are based on similar kind of assumptions as used for Glauber approach used in heavy ion collisions.
Initial transverse shape of the nuclei as described by Glauber model for heavy-ion collisions depends on Wood-Saxon distribution, which is a twoparameter (half-density radius (R) and diffusivity (a)) Fermi-like distributions (2pF) extracted from fits to elastic lepton-nucleus data [21,22], which describes the multi-nucleon interactions occurring in the overlap region between the colliding nuclei via a Glauber eikonal approach [23]. Whereas, in the Monte Carlo Glauber (MCG) models [24][25][26][27][28][29], event-by-event sampling of individual nucleons are done from Wood-Saxon distribution and averaging over multiple events are used to calculate properties related to collisions. Presently, available partonic Glauber model for p+p collisions does not consider full anisotropic density profile of protons, though radial homogeneity is assumed.
In this article, we present the results of Glauberlike model calculations for N coll (b), N part (b) due to the quark and gluon based proton density profile, which is a realistic picture obtained by results of deep inelastic scattering that reveals the structure of proton [31], and used it to obtain charged particle multiplicity distribution in p + p collisions at √ s = 7 TeV. Calculated multiplicity distribution is contrasted with ALICE data, a relation of impact parameter with multiplicity is calculated and multiplicity distribution of eccentricity and flow harmonics is estimated for p+p collisions. In order to understand the possibility of medium formation in highmultiplicity p + p collisions, we have estimated nuclear modification-like factor, R pp , considering low multiplicity yields as the base.
The paper is organised as follows. In the next section we discuss the formalism that is used in this work. In section III, we present the results and section IV is devoted to summary and discussions.

II. GLAUBER FORMALISM
In the literatures, the density profiles like hard sphere, and 2pF functions are used traditionally to formulate Glauber model for heavy ion and even for protons [29]. All these profiles can also be extended to proton model by considering radially symmetric parton density. In fact, in the case of proton, several density profiles have been considered to estimate the initial conditions, most of them assume azimuthally symmetric density profile, those are mainly different in phenomenological paramatrization of radial variations [30]. But the standard model postulates that a proton consists of three effective quarks (constituent) and gluons within it. Thus distribution of such configuration is less likely to be radially symmetric, because we expect individual peaks in wave function in the quarks position inside a proton indicating its presence. The necessary condition is, however, that the wave function of each effective quarks and gluons should decay rapidly around boundary of a proton (within RMS area). In this regard we find only one previous work [31] to consider azimuthally asymmetric and homogeneous density distribution of proton, which is motivated by the shape of deep inelastic scattering structure functions, pointing out that multiplicity distribution produced by different model can be used to discriminate them-which can better reproduce experimental results.
In this study, we have used a model with fluctuating proton orientation and it has three effective quarks and gluonic flux tubes connecting them. The densities of quarks (ρ q ) and gluons (ρ g ) are taken as Gaussian type assuming spherically symmetric distribution of quark densities from their respective centers and cylindrically symmetric gluon densities about the line joining two adjacent quarks as where, r q is the radius of a quark, r s and r l are the radius and the length of the gluon tube. The density function under study here was taken to be [31], where, R[θ, φ] transforms vector (0,0,1) into (cosφsinθ, sinφcosθ, cosθ) and r k = r k (cosφ k sinθ k , sinφ k cosθ k , cosθ k ) (where, k = 1,2 and 3) is the position vector of k th effective quark. N g is the number of partons inside a proton. The free parameter κ allows to control the percentage of gluon body content and here it is taken to be 0.5 [31].

A. Calculation of Thickness function and Overlap function
The collision plane is taken to be in x-y hence dependence along z axis is integrated out as follows The calculated thickness function for the ρ G−f is and and The overlap function T pp (b) for projectile proton (A) and target proton (B) is defined as Here T pp is sum of 4-components namely quarkquark, quark-gluon, gluon-quark, gluon-gluon. Primed (unprimed) indices indicate variables corresponding to B (A). In the following, we provide overlap function for all the possible combination of partons.

A.1. The Quark-Quark term
The overlap function for the interaction of two quarks:

A.2. The Gluon-Gluon term
The overlap function for the interaction of two gluon tubes: where, The overlap function for the interaction of a quark and a gluon tube: where, A.4. The Gluon-Quark term The overlap function for the interaction of a gluon tube and a quark: where, Together total overlap function is sum of four terms given by Eq. 12, 13, 18, 23.

B. Calculation of N coll and Npart
We define the number of binary collisions (N coll ) of partons in a p + p collision at a given impact parameter (b) as follows: In line with the previous studies [31,32], we fix σ gg = 4.3 ± 0.6 mb [33] with N g = 10 partons, so as to reproduce the experimental value of inelastic cross section, σ pp = 60 mb for p + p collision at √ s = 7 TeV. This accounts for the only non-trivial dependence of the Glauber calculation on the beam energy √ s. Previous studies [31,33] have assumed linear scaling of charged hadron (N ch ) multiplicity with N coll only. In contrast to this assumption, we have considered dependence of N ch on number of participants partons (N part ) and N coll . Further, relationship between N part and N coll is considered non-linear as that of the heavy ion collisions assuming three dimensional shape. Thus, number of participating partons at impact parameter 'b' is given as where x is a parameter.
By considering, f , as a fraction of charged hadron multiplicity produced from binary collisions, we have two component model for estimation of number of charged particles given as where, n pp is a constant of proportionality and f is a free parameter.

III. RESULTS
Assuming initial position vectors of three quarks to be vertices of equilateral triangle in xy-plane as r 1 = ( d 4 , where, d is the free parameter of the model which ensures that the length of the gluon tubes connecting quarks are fixed i.e., (|r 1 | 2 = |r 2 | 2 = |r 3 | 2 ) = d 2 4 ). Now, in order to account for all possible configurations, position vectors of quarks are parameterised by varying azimuthal and polar angles. Generalised configuration considering tilt by ψ along x-axis and rotation by angle α are sin ψ sin α) and considering tilt by γ along y-axis and rotation by angle β are r 1 = ( d 2 cos( π 3 + γ) cos β, d 2 sin( π 3 + γ), d 2 cos( π 3 + γ) sin β), In the above configurations, ψ and γ (0, 2π 3 ), α (0,π) and β (0,2π). In our present study, we have taken, x, in Eq. 30 to be 0.75 as N coll scales as A 4/3 for similar target and projectile nuclei with mass numbers A for heavy ion collisions and are spherical in shape. In our work, this consideration of x = 0.75 holds good because when the plane formed by connecting centres of each quark is randomly rotated as part of Monte Carlo simulation for accounting all possible configurations of collision geometry, the overall angular space is exhausted thus making collisions geometry to be closely spherical overlap with preserving contributions from each of the different configurations hence the factor 0.75 is taken so that it accounts for general spherical overlap in heavyion collisions. We have also chosen, RMS radius of proton and quark as 1 fm and 0.25 fm, respectively.

A. Number of binary collisions and participants as a function of impact parameter
We have used Eqs. 29 and 30, to estimate N coll and N part . Fig.1 shows the mean value of N coll (upper curve) and N part (lower curve) as a function of impact parameter (b). Towards higher values of b, the difference between the two curves effectively vanishes. Similar trends were observed for Au+Au and Cu+Cu collisions at √ s N N = 200 TeV [19].

B. Charged particle multiplicity estimation
Two-component models have been used in heavyion phenomenology for long time to estimate the charged-particle multiplicity [34,35]. The inelastic cross section, σ inel N N , which depends on collision energy, is used as input for the MC Glauber model. In our current study, we have used similar approach for p + p collisions as well, where nucleons are replaced by partons (quarks and gluons) and σ inel N N by σ inel gg . The model provides N part and N coll , for an event with a given impact parameter and collision energy which is discussed in the previous section. As in heavy-ion collisions, the concept of "ancestors" (independently emitting sources of particles) has been introduced for a given value of N part and N coll . The number of ancestors can be parametrized by a twocomponent model given by [34,35], The two-component model divides the partonparton collisions into sof t and hard interactions: the multiplicity of particles produced by soft interaction is proportional to N part and hard interaction is proportional to N coll . As negative binomial distribution (NBD) able to well reproduce the chargedparticle distribution in p + p collisions [36], we use two-parameter NBD to calculate the probability of producing n particles per ancestor; P (n;n, k) = Γ(n + k) Γ(k)Γ(n + 1) n k +n n k k +n k , (33) wheren is the average multiplicity and k characterizes the width of the distribution. By the use of different combination of f (Eq. 32),n and k ( Eq. 33) we have repeated the process of obtaining the multiplicity distribution for a large sample of events, until our model simulate the experimental multiplicity distribution. We have also calculated the ratio of N ch obtained from our model to that of experimental value and is represented in Fig. 2 for p + p collisions at √ s = 7 TeV. The best agreement for N ch distribution obtained by our model with experimental data is found for f = 0.85,n = 8 and k = 0.13. From Fig. 2, it can be seen that our model well describes the data in the mid multiplicity region ( 15 < N ch < 90 ), with 5-10 % discrepancy. However, towards the low and high multiplicity it is unable to reproduce the experimental measurement. The inability of the model to explain the extreme low and high multiplicity region might be due to lack of statistics.

C. Centrality estimation
The centrality is usually expressed as a percentage of the total interaction cross section, σ [38]. Impact parameter distribution is taken as input to our current model. So, the centrality percentile of a p + p collision with b is defined by integrating the impact parameter distribution as,  1 , b 2 ) and so on. For the current analysis, a Gaussian distribution with mean 1 and standard deviation of 0.32 has been used as input impact parameter distribution, which is shown in Fig. 3, so that the distribution function vanishes beyond the proton radius (≈ 2 fm).
We have also tested different forms of impact parameter distributions, but Gaussian distribution is found to be a suitable choice to describe the chargedparticle multiplicity distribution. Once, we get the ranges of the impact parameter corresponding to each centrality, we have projected it to N ch , N part b (fm) and N coll to calculate N ch , N part and N coll corresponding to each b-ranges. Fig. 4 represents the multiplicity distribution for each percentile bin. Table I shows the value of N ch , N part and N coll , calculated using our model along with N ch value of ALICE for pp collisions at √ s = 7 TeV. It can be clearly seen that the calculated dN ch /dη is well consistent with experimental value, except for the high and low multiplicity regions. This is because of the artefact of incapability of our model to describe the charged-particle distribution in that region (Fig. 2). However, it is to be noted that the input σ inel gg = 0.43 ± 0.06 f m 2 contain 14 % uncertainty and the same amount of uncertainty (14 %) is associated with each dN ch /dη . From our model, we found dN ch /dη = 7.47 for minimumbias (0-100 %) collisions, which is little higher from the experimental value, dN ch /dη = 6.01 ± 0.01 +0. 20 −0.12 [39]. This discrepancy needs to be understood.
D. The ratio, Rpp for high to low multiplicity events In order to understand the possibility of formation of a medium in high-multiplicity events in p + p collisions, we define a variable as: which is similar to the nuclear modification factor R AA in heavy-ion collisions. Here, d 2 N/dηdp T | HM , d 2 N/dηdp T | LM , N LM coll ( N HM coll ) are charged particle yields in high-multiplicity, low-multiplicity p+p collisions at √ s = 7 TeV [40], mean number of binary collisions in low (high) multiplicity p+p events. Upper panel of Fig. 5 shows the transverse momentum spectra of charged particle in high-multiplicity (VOM I), second high multiplicity (VOM II) and low multiplicity (VOM X) events obtained from Ref. [40]. And lower panel shows the R pp defined in Eq. 35. For such definition of R pp , it is observed for all charged particles for p T < 1 GeV/c, value of R pp < 1 and for p T > 1 GeV/c , it is greater than 1. However, it tends to reduce at very high p T . And for p T >1 GeV, the value of the factor is higher for higher multiplicities. Fig. 6 shows results of R pp for identified particles, pion (π + + π − ), Kaon (K + + K − ), proton (p +p) for p + p collisions at √ s = 7 TeV [40]. It is found that R pp < 1 for proton for p T < 1 GeV which is same as observed in case of charged particles. However for pion and kaon R pp < 1 for p T < 0.8 GeV. It is also observed that for p T < 1.9 GeV, this identified particles have almost value of R pp and for p T > 1.9 GeV, the value is almost same for pion and kaon but the value for the proton is much more and increases with p T sharply upto p T = 5 GeV, and then saturates within uncertainties. But for pion and kaon, the factor increases monotonically with decreasing slope from p T > 1.9 GeV, where the trend splits for proton and other two hadrons.
It is reported that [41], proton shows distinct behaviour in this regard than other hadrons produced in p-Pb collisions. Also for p-Pb collisions, it is reported that the factor, R pP b > 1, for all charged particle for p T > 2.5 GeV [41,42]. For p-Pb, R pP b saturates to unity for p T > 2 GeV, and here it is found that for p + p, R pp shows almost similar trend  [40] for VOM multiplicity classes, viz., highest (HM), second highest (second HM) and lowest multiplicity (LM) class. Lower Panel: Rpp obtained from the ratio of differential yield at high-multiplicity and second high multiplicity classes with low multiplicity class scaled by N coll . but with larger value of the factor with saturationlike behaviour starting after p T = 2 GeV. We note that the R pp values above unity for p T > 1 GeV may be qualitatively similar to other observed enhancements due to the Cronin effect and radial flow in pA and dA systems [43,44], as conjectured for similar behaviour of R pP b [42], where the moderate excess at high p T is suggestive of anti-shadowing effects in the nuclear parton distribution function [45].

E. Estimation of Elliptic-flow
For long time, p + p collisions were considered as the baseline measurements for determination of deconfined state of matter i.e., QGP. Recent observation of p + p collisions at LHC energies hints toward collective effect, thus, it becomes evident that the earlier view needs upgradation. In this regard, we have also calculated eccentricity ( ) using the present approach. Asymmetry ratio between semi-axis dimensions of the overlap region weighted by N coll at a particular b can be used to obtain as: where, n coll (x, y, b) = σ gg T a (x − b 2 , y)T b (x + b 2 , y) represents impact plane binary collision density. We have calculated (b) by using Eq. 36 by considering sum of 4-components namely quark-quark, quarkgluon, gluon-quark and gluon-gluon. Fig. 7 shows the eccentricity for p + p collision at √ s = 7 TeV obtained using Eq. 36 and it is observed to increase with b and seems to saturates towards larger b.
Using , we have obtained elliptic flow (v 2 ) as a function of b by considering v 2 scaling i.e., v 2 = Ω , where Ω = 0.3 ± 0.02 [33]. By geometry, v 2 (b) will follow the general trend of (b). It is found that the overlap of two hard spheres with infinitely sharp edges yields artificially large eccentricities [49].
In Fig. 8, we have compared our estimation of variation of v 2 with charged particle multiplicity for p + p collision at √ s =7 TeV with experimental result at √ s =13 TeV [50]. This is due to the fact that, when reporting this work, we do not have the data for collisions at √ s =13 TeV to constrain our model, for that we have only data for collisions at √ s =7 TeV. That does not prevent us from the comparison, since in Ref. [51], it is reported that value of v 2 for collisions at √ s =2.76 TeV and √ s =13 TeV are almost the same when measured for different transverse momentum, indicating v 2 is independent of collision energy. It is observed that for N ch 8, our estimation of v 2 with linear response to initial geometry reproduces the value obtained from experiment within the error bars. However, for lower multiplicities, our estimation with linear response to initial eccentricity falls short to that obtained from experimental data. This may be due to effects other than collective linear response or final state effects. Though, the charged particle multiplicity variation of v 2 for p + p collisions at √ s = 7 TeV is not available, the elliptic flow coefficient extracted from the CMS Collaboration data at √ s =7 TeV is 0.04 − 0.08 [52] and our estimation of v 2 falls within this range. We also note that this model gives v 2 similar to that of the IP-Glasma model as presented in Ref [53] for low multiplicity region (< 8).

IV. SUMMARY AND DISCUSSION
In this work, we have investigated predictions of Glauber model for the initial condition for p + p collisions, which considers anisotropic and inhomogeneous proton density profile, by contrasting with experimental data. This model for density profile is inspired by the structure function obtained from deep inelastic scattering. Instead of distributing the positions of valence quarks randomly by keeping center of mass intact, we have taken random orientations generated by random rotation around three spatial axes, where center of three quarks form a plane and connecting gluon tubes always remain fixed in length. This prevents the overlap of two valence quarks in space and possible placement of a quark out of proton radius, where these two can happen for the first kind of randomization with only center of mass being fixed [31], and to condition which may bring extra complication in the randomization process for not allowing it, generating "spooky" correlations. However, the present approach, apart from avoiding such complications, will give better handle for future investigations.
With all these considerations, we have studied multiplicity distribution, to obtain the impact parameter to multiplicity relation, multiplicity dependence of initial eccentricity and azimuthal flow harmonics (v 2 ). It is found that this model can well reproduce multiplicity distribution produced in experimental p + p events at ALICE, with the free parameter f = 0.85. With properly constraining our model with experimental data and calibrating the range of b with multiplicity percentile, we have used the estimated N coll to obtain nuclear modificationlike factor (R pp ) for p + p collisions. It is found that the defined factor < 1 for p T < 1 GeV, and beyond this, the factor > 1. Moreover, it tends to reduce at very high p T , and for p T > 1 GeV, the value of the factor is higher for higher multiplicities. We have also studied R pp for identified particles for p + p collisions at √ s = 7 TeV, and found that the trend for R pp is similar to that observed in p-Pb system but with increased value. This behaviour at higher p T may be due to non-collective flow effects, which needs further investigation.
The non-availability of results from experiments which shows the variation of eccentricity and v 2 with multiplicity at √ s = 7 TeV prevents us from comparing our estimation with experimental data at √ s = 7 TeV. However, we have compared our result of v 2 with that of p + p collisions at √ s = 13 TeV, as it is observed that the collision energy dependence of v 2 is weak. We found result of v 2 obtained from present approach to be in agreement with IP-Glasma model for lower multiplicity region. Also, it is found that the values of v 2 obtained from present model for N ch 8 are very close to that of experimental data for √ s = 13 TeV. Such good agreement with data for multiplicity distribution and results of v 2 will make it a reliable initial condition to investigate fluid nature of system produced in p + p collisions. Kubiczek for discussion on geometrical models for p + p collisions.
VI. APPENDIX Table  I shows geometric properties ( b , N ch , N part , N Coll ) of p + p collisions for different multiplicity classes using Glauber Monte Carlo calculation along with a Negative binomial distribution fit to charged particle multiplicity distribution at √ s = 7 TeV for the ALICE experiment at the LHC.