Investigating high energy proton proton collisions with a multi-phase transport model approach based on PYTHIA8 initial conditions

The striking resemblance of high multiplicity proton-proton (pp) collisions at the LHC to heavy ion collisions challenges our conventional wisdom on the formation of the Quark-Gluon Plasma (QGP). A consistent explanation of the collectivity phenomena in pp will help us to understand the mechanism that leads to the QGP-like signals in small systems. In this study, we introduce a transport model approach connecting the initial conditions provided by PYTHIA8 with subsequent AMPT rescatterings to study the collective behavior in high energy pp collisions. The multiplicity dependence of light hadron productions from this model is in reasonable agreement with the pp $\sqrt{s}=13$ TeV experimental data. It is found in the comparisons that both the partonic and hadronic final state interactions are important for the generation of the radial flow feature of the pp transverse momentum spectra. The study also shows that the long range two particle azimuthal correlation in high multiplicity pp events is sensitive to the proton sub-nucleon spatial fluctuations.


Introduction
Collective phenomena have been regarded as important signatures indicating the creation of a new state for nuclear matter, Quark-Gluon Plasma (QGP), in relativistic heavy ion (AA) collisions [1,2,3]. In recent years, a flood of striking collectivity like features have also been observed in high multiplicity proton-proton (pp) and protonnucleus (pA) collisions at the Relativistic Heavy-Ion Collider and the Large Hadron Collider (LHC) [4,5]. These observations were mostly unanticipated for small systems like pp which were expected to be QGP free references for the study of the de-confined quark gluon matter produced in AA collisions. An early example like this comes from the observed ridge structure of the azimuthal correlations of two charged hadrons with a large pseudo-rapidity gap [6,7,8]. Later discoveries on the enhancement of multistrange particles [9,10,11], non-vanishing elliptic flow coa e-mail: zhengliang@mail.cug.edu.cn b e-mail: linz@ecu.edu c e-mail: shouqiye@fudan.edu.cn d e-mail: zbyin@mail.ccnu.edu.cn efficients [8,12,13,14], mass dependence of the hadron p T spectra [15,16,17,18] and characteristic variations of the p T dependent baryon to meson ratios [17,19,16,20] further consolidate the resemblance between the small system and AA events.
The observation of the surprising collectivity in pp collisions challenges our conventional wisdom on the formation of the QGP medium and invites a lot of theoretical studies to understand its origin [5,21,22,23]. It is conceived that a QGP droplet can be partly formed in high multiplicity pp events and experience the same hydrodynamic expansion as in AA collisions [24,25,26,27,28,29]. Similar ideas are also implemented with a corecorona picture in the EPOS model [30,31,32]. Alternative approaches focusing on the initial state gluon correlations based on the color glass condensate are also found to be successful in describing the qualitative features of the seemingly collective behavior [33,34,35,36].
Another group of works relies on extensions to the traditional string fragmentation model in the PYTHIA event generator [37,38] by considering interactions between multiple string objects when they overlap in the transverse arXiv:2104.05998v2 [hep-ph] 27 Aug 2021 space. The color reconnection model revises the string structure created in the events with multiple parton interactions (MPI) and is useful in dealing with the flavor dependent observables [39,40]. The color rope model, initially applied to pp collisions in DIPSY [41], now a part of the standard PYTHIA program, is also capable to reasonably describe the multiplicity dependent particle compositions [42]. However, dynamical correlations, like the near side ridge and elliptic flow observables, are generally difficult to reproduce within the traditional string fragmentation framework. The recent development of the string shoving mechanism [43,44,45] is expected to improve the description of these measurements by considering the repelling of string pieces overlapped in space-time.
On the other hand, the string melting version of a multi-phase transport model (AMPT) [46], which includes both partonic and hadronic final state rescattering effects, has been extensively applied to interpret the collective flow both in AA and in small systems [47,48,49,50,51,52]. With the event by event geometry fluctuations, this method provides a generic way to understand the importance of the parton degrees of freedom in collectivity observables and shows the non-equilibrium effect such as the parton escape mechanism can be quite significant for the development of collective flow especially in small systems [22,23,53]. Whilst being successful in depicting the azimuthal anisotropic flow, the AMPT model can not easily reproduce the multiplicity dependence of transverse momentum spectra until the latest improvement [54]. In this work, we build a transport model framework based on the initial conditions provided by PYTHIA8 linked to the final state interactions and quark coalescence model [55] within AMPT, in an effort to simultaneously describe the p T spectra and multi-particle correlations in pp collisions at √ s = 13 TeV. One can exploit the interplay of parton and hadron final state interactions by considering the microscopic dynamical process of the evolving system. This study may pave the way for further investigations on the azimuthal anisotropy measurements in small systems and help us to fathom their origins from the perspective of a transport model. The rest of this paper is organized as follows: we elaborate the transport model approach for pp collisions and its extensions with the proton transverse spatial geometries in Sec. 2. Comparisons of the model calculations to the multiplicity dependent pp data from the LHC are made in Sec. 3. Finally, we summarize our major conclusions and discuss the potential applications of this framework in Sec. 4.

Model setup
Measurements in small collision systems have shown anisotropic flows similar to those observed in heavy ion collisions in which hydrodynamic descriptions are found to work well [56]. This implies that a deconfined quark gluon matter may have been created in these small systems. In the PYTHIA framework, the partonic system including jets and perturbative radiations is kept in the strings even in the high energy density region, which might underestimate the collective partonic effects. To account for the effects of the parton degrees of freedom, we adopt the string melting mechanism from AMPT, converting excited strings generated by PYTHIA to their constituent partons, and couple the resultant parton system to further final state rescatterings based on the kinetic transport model approach.
A high energy pp event in this approach is described by the same four ingredients as the AMPT model: initial conditions, partonic rescatterings, hadronization and hadronic rescatterings. Unlike the full AMPT model which uses HIJING [57] to obtain initial conditions, PYTHIA8 is utilized to generate the initial conditions for the subsequent evolution stages. The PYTHIA/Angantyr model [58] implemented in PYTHIA8 can provide an impact parameter dependent string system within the MPI framework suitable for further treatment considering its space-time structure. We consider the primary hadrons from the excited strings as the intermediate step to model the string melting process. The parton system is obtained by converting the primary hadrons to their valence quarks [59]. A formation time is given to each primary hadron including its valence quarks, with E H and m T,H being the energy and transverse mass of its mother hadron. Partons do not have any interactions until after their formation time. The primary hadrons can be regarded as the parton sources, which give the parton initial space-time configuration after the propagation of the primary hadrons to their formation times. The microscopic dynamic processes of the parton evolution are implemented with the Zhang's Parton Cascade (ZPC) model [60] by solving the Boltzmann equation of two-body elastic scattering process through the parton cascade method. The parton rescattering cross section in ZPC determines the evolution feature of the deconfined parton medium. For the partons without further scatterings in ZPC, they will become hadrons through a spatial quark coalescence model. The hadron species are determined by the coalesced (anti)quark flavor combination with an overall parameter r BM = 0.53 dictating whether a meson or a baryon is formed [55]. Further hadronic scatterings are treated with an extended relativistic transport model (ART) [46,61] with both elastic and inelastic scatterings.
The key model parameters used in this approach can be divided into two categories, the parameters involved in the initial conditions and those relevant for the final state rescatterings. As the initial conditions are generated by PYTHIA8 in this framework, we take the Monash tune of PYTHIA8 parameters in the simulation. The key final state interaction parameter is the parton rescattering cross section in ZPC. Its value is determined by systematic comparisons with the elliptic flow data. The full AMPT model with the new quark coalescence uses 1.5 mb for large systems [55], while in this work we provide results with 0 mb, 0.2 mb and 3 mb in the comparison to demonstrate the parton rescattering effects based on our current approach. If we turn off all the final state interactions, the result is expected to be very similar to PYTHIA8 Monash calculations.
Subnucleonic geometry has been speculated to be important in the space-time development of the evolving nuclear matter produced in small collision systems. In our transport model approach, we take the sub-nucleon geometry into consideration when sampling the initial transverse positions of the parton sources before converted to constituent quarks. The sub-nucleon geometry may reveal itself together with strong spatial fluctuations in small system collective flow [56,62,63,64,65]. A precise constraint on the underlying parton spatial distribution of the nucleon might require the data from future deep inelastic scattering experiments [66]. We perform our study with two initial geometry models representing our current understanding of the proton inner structure.
Different matter distributions for proton have been extensively compared in [67] and considered to be strongly connected to the eccentricity of the evolving system produced in pp events. It is common to parameterize the proton spatial density ρ(r) with the following different distributions: hard sphere, exponential, Fermi, Gaussian and double-Gaussian functions. We provide results for the exponential matter distribution based on the proton charge form factor where R = 0.2 fm. The transverse coordinates of the produced string objects are determined according to the overlap function of a pp collision at an impact parameter b with T (x, y, b) = ρ p,1 (x − b/2, y, z)ρ p,2 (x + b/2, y, z)dz assuming the two colliding protons are moving along the z axis. A schematic illustration of the spatial distribution of the produced objects based on the overlap function T (x, y, b) weighting method is provided in Fig. 1 (a). The previous method gives a smooth initial transverse spatial condition but may underestimate the fluctuations in the proton matter distribution. We also introduce a Glauber Monte Carlo type method considering the event by event sub-nucleon level fluctuations based on the constituent quark picture [68,69,70,71]. By assuming that a proton can carry three constituent quarks, the interaction of two protons can be extended to the participant quark geometries within the Glauber model framework. The constituent coordinates in a proton can be sampled according to the proton form factor as shown in Eq. 2. Two constituents from the different incoming protons thus may collide if their transverse separation is less than D = σ cc /π, with σ cc being the constituent scattering cross section. The transverse coordinates of the parton sources are randomly assigned to the binary collision center of each interacted constituent pair, as shown in Fig. 1 (b). The results shown in this work are obtained with σ cc = 25.2 mb, which reproduces the inelastic pp cross section at √ s = 13 TeV [69]. It must be emphasized that the constituent quark picture used in our current approach is different from an initial state quark participant Glauber Monte Carlo model [71]. We first take the interacting nucleon collision geometries provided by PYTHIA/Angantyr for pp and then assign the produced particle transverse positions using the colliding quark spatial conditions as an afterburner to account for the event by event sub-nucleon fluctuations.
As a comparison, the initial transverse spatial distributions of the parton sources produced in AMPT is illustrated in Fig. 1(c). The parton system produced in AMPT based on the HIJING initial conditions can be divided into two types. The first type coming from the wounded projectile or target proton is attached to the position of the incoming projectile (b/2, 0) or target (−b/2, 0) beam in transverse plane, respectively. The second type from the possible minijet production process is put at the midpoint (0,0) of two beam locations. Unlike the approach in our current work modeling the parton spatial distributions based on the proton geometry in a dynamical way, the full AMPT model is placing the initial parton objects at three fixed locations largely dependent on the size of the impact parameter.
In kinetic transport model methods, the final state particle momentum anisotropy is converted from the initial spatial structures via the final state rescattering process. The initial spatial structure of an event is tightly related to the impact parameter distribution and the parton number created in the evolving system. We present the impact parameter distribution taken from PYTHIA/Angantyr for pp 13 TeV collisions in Fig. 2(a). It can be observed in this figure that the maximum impact parameter may rise up to 5 fm. The parton number distribution created after string melting process over the full phase space is shown in Fig. 2(b). The number of created partons in pp collisions at the LHC energy can be very large.
We present the initial parton spatial distributions in the transverse plane right after string melting in Fig. 3 for different proton geometry assumptions. Only partons with formation time t p < 5 fm/c are displayed in this figure. The solid dots represent the initial parton positions at each parton's individual formation time, while the thin arrows indicate their transverse velocity vectors (with the length proportional to the velocity magnitude β ⊥ = p T /E). The thick red arrow gives the direction of the second order event plane obtained from the positions of the shown partons. The spatial eccentricities and event plane directions are obtained following the calculations made for participating nucleons in Ref. [72]. Partons that will experience rescatterings in ZPC are shown in red and the light green partons have no interactions in parton cascade. The two black circles with radius R = 0.8 fm demonstrate the transverse view of the proton beams located at (−b/2, 0) and (b/2, 0). The first row represents results obtained with initial conditions based on the T (x, y, b) weighting method and the second row includes the subnucleon fluctuations within the constituent quark picture. We find that the transverse spatial structure of the produced initial parton system largely depends on the implemented proton geometry. With T (x, y, b) weights, the partons are mostly distributed around the center of the overlap region for two proton beams at b = 0.1 fm as  shown in Fig. 3 (a). In peripheral collisions (b = 0.6 fm), the parton system becomes stretched along the impact parameter direction in Fig. 3 (b). On the other hand, it is observed in Fig. 3 (c) and (d) that multiple hotspots or localized high density regions can be produced with the constituent quark picture around the binary collision centers of each interacted quark pair shown by the yellow stars. Figure. 4 shows the average initial eccentricity distributions obtained from the parton spatial positions right after string melting versus the parton number and impact parameter. Similar to Fig. 3, only partons with t p < 5 fm/c are considered for the estimation of eccentricity while N parton in Fig. 4(a) represents the number of partons in the full phase space melted from the initial primary hadrons in that event. Eccentricities are calculated with the initial positions at each parton's formation time [50]. For the overlap function weighting method, the eccentricity is largely driven by the geometric shape of the transverse overlap area. The eccentricity drops significantly from peripheral collisions to central collisions. However, the quark constituent picture that has sub-nucleon fluctuations generates larger eccentricities in central collision events in contrast to the overlap function weighting method. The impact parameter dependence of the eccentricities from the quark constituent model is quite weak. We also present the eccentricity distribution from the full AMPT model in Fig. 4 with the blue dotted line. Its value is found to be significantly larger than that from our current approach except for the events with very small b. The origin of the eccentricity in AMPT and how it affects the final state collectivity will be discussed at the end of Sec. 3.
With the overlap function weighting method, time evolutions of the energy density and the number density of the partons are shown in Fig. 5. In this comparison, the parton densities are obtained with all the active partons which are formed at t p < t and not yet hadronized at time t. The density distributions of the parton system in the central cell of pp collisions at √ s = 13 TeV with b = 0 fm using different final state parton rescattering cross sections are presented. The density of partons decreases at large t due to the expansion of the system. The solid line represents the transverse energy density ε T and the dashed line shows the parton number density n. The central cell is defined with a cross sectional area with the transverse radius of 0.5 fm and a longitudinal dimension covering −0.5t to 0.5t. The energy density of the system is found to be above the critical energy density (∼ 1 GeV/fm 3 ) within t < 1 fm/c, indicating the possibility of producing a short lived QGP droplet in high energy pp systems. By increasing the parton rescattering cross section in ZPC from 0.2 mb to 3 mb, the energy density of the system be- comes larger at any given time during the entire evolution stage. This observation implies that the deconfined matter in pp collisions with stronger final state interactions can have a longer lifetime. It is also worthwhile to note that the typical formation time of the partons is around 0.1-0.5 fm/c in the string melting scheme, different from the string fragmentation time at about 1 fm/c. As is emphasized in Ref. [46,50], the general time scale in the string fragmentation picture may not apply to the initial partons as string melting mechanism models parton not hadron productions from the string energy field.

Results
The final state particle productions are compared to the pp 13 TeV experimental data after all the secondary interactions and resonance decays are finished in the transport model. It is found in our study that different proton geometries for the initial condition have a negligible impact on the inclusive particle productions. Initial state conditions mainly reveal themselves in the correlation measurements. Therefore, we will show the effects of final state interactions on particle yields and spectra using the constituent quark initial condition in Sec different initial conditions on long range two particle correlations in Sec. 3.2. The results from the public AMPT code (v2.26t9) [73] are also presented in the comparisons with Lund fragmentation parameters taken as a = 0.5, b = 0.9 GeV −2 and the parton rescattering cross section set to 1.5 mb.

Particle spectra and final state interactions
The charged hadron pseudo-rapidity distributions are investigated in Fig. 6  actions are shown in dashed and dotted lines, respectively. A quite large parton rescattering cross section 3 mb has been used in this comparison to explore the maximized parton evolution effects. Calculations are compared to the minbias measurements from ALICE inelastic events (black markers) and INEL>0 events with at least one charged track in mid-rapidity (red markers) [74]. It is interesting to observe that, unlike AA collisions, final state interactions do not change the charged multiplicity density too dramatically in pp. As parton rescatterings generally decrease the parton transverse momentum, it is observed that the mid-rapidity particles are shifted to forward and backward η regions by comparing the solid and dashed lines. On the other hand, the transition from solid to dot- ted lines is an indication of the radial flow effect induced by the hadronic rescatterings, which causes the excess in the charged particle density at η ∼ 0. Stronger radial expansion arised from hadron rescatterings increases the average p T of charged hadrons and pushes more particles to the mid-rapidity region. It is also found that the long dahsed line, with combined parton and hadron final state effects and a moderate parton rescattering cross section 0.2 mb, reproduces the INEL>0 charged particle density data reasonably, although the inelastic event charged particle yields are overestimated.
In Fig. 7, we compare the transverse momentum spectra for π ± , K ± and p(p) in our model to the ALICE pp data at √ s = 13 TeV [16] in different centralities. The percentile definition of the events are made with the charged track number in the pseudo-rapidity regions −3.7 < η < −1.7 and 2.8 < η < 5.1 following the event classification scheme used at the ALICE experiment. Ten event classes (Class I ∼ Class X from central to peripheral) are obtained with the same analysis procedure in [16]. The results for the most peripheral event class 64.5-100% (Class X) and the most central event class 0-0.92% (Class I) can be found in Fig. 7 (a) and (b), respectively. Similar to what has been done in Fig. 6, we only show the comparisons of final state effects based on the constituent quark proton matter distribution. Considering that the string melting framework generally works better for low to intermediate p T region, we focus on the p T range from 0 to 3 GeV/c in the spectra comparisons. Results without final state interactions are presented in solid lines. Short dashed and dotted lines represent the calculations with only parton rescatterings (3 mb in ZPC, ART turned off) and only final hadron interactions (0 mb in ZPC, ART turned on), respectively. The long dashed lines represent the calculations that include both parton rescatterings and final state hadron interactions. π ALICE ALICE K ALICE p Fig. 9. Average transverse momentum of π (black), K (blue) and proton (red) at mid-rapidity with 0 < pT < 3 GeV/c versus the charged hadron multiplicity density. Final state effects are shown in (a). Comparisons to AMPT are shown in (b). Central values of ALICE data are obtained with a refit to the data from [16]. In the peripheral event class with low multiplicities shown in Fig. 7 (a), the final state interaction effects are expected to be very weak. Differences between various final state interaction combinations are indeed found to be small. In central pp collisions shown in Fig. 7 (b), a significant stiffening can be observed for the p T spectra with only hadron rescatterings effects (dotted line). The scenario of pure hadronic rescatterings has a strong radial expansion, boosting the hadrons to higher p T with a mass ordering feature. Protons receive stronger change of p T compared to the lighter K ± and π ± . On the other hand, parton rescatterings are expected to decrease the hadron yield at larger p T . When parton and hadron interactions are considered together, the strong radial expansion observed in the pure hadronic rescattering scenario is no longer there (long dashed line) because the parton stage delays hadronic rescatterings to lower densities. Underestimations in high p T region exist for π ± and K ± but not for p(p) since the proton transverse momentum suppres-sion due to parton rescatterings is compensated by the stronger radial flow from secondary hadronic interactions.
We also compare our current results on the multiplicity dependent transverse momentum spectra to the calculations from the public AMPT code in Fig. 8. In peripheral pp collisions, only the proton yield is slightly overestimated in AMPT and the meson p T spectra are generally well reproduced. In central pp events, the AMPT model predicts softer transverse momentum spectra especially for protons. Compared to the full AMPT model, our current approach improves the description of the identified particle p T spectra in high multiplicity events. Figure 9 (a) shows the final state effects on the multiplicity dependent average transverse momentum distributions for π, K and p(p). Only hadrons with 0 < p T < 3 GeV/c at mid-rapidity are considered in this comparison. The ALICE p T data are obtained by refitting the measured p T spectra with the Lévy-Tsallis parameterization restricted to the corresponding p T range while the uncertainties are taken from the experimental data with full p T range [16] directly. Similar to the observation in Fig. 7, the scenario of pure hadronic rescatterings shown by dotted lines leads to a strong enhancement in the average p T . The magnitude of the enhancement is mass dependent and more pronounced for protons. The scenario of pure parton rescatterings (short dahsed lines) suppresses p T for pion and kaon but slightly enlarges proton p T . Figure. 9 (b) shows a comparison of our current model calculations with all final state interactions to the full AMPT results together with the experimental data of p T . By including both the parton and hadron evolutions, results from our current model calculations reproduce the trend of the data within the selected p T range, while overestimations are found for π, p(p) average p T and kaon p T is slightly smaller than data. The full AMPT calculations reasonably describe the pion p T but generally deliver very weak multiplicity dependence for kaon and proton. By decreasing the parton rescattering cross section in the full AMPT from 1.5 mb to 0.2 mb, it is found that the average p T only slightly changes and the flat multiplicity dependence persists. The improvement compared to the full AMPT model on the multiplicity dependence of kaon and proton p T presumably comes from the MPI mechanism [75,76] implemented with the PYTHIA initial conditions incorporated to our transport model framework.
We present a group of comparisons of the final state interactions on the baryon to meson ratio p/π in Fig. 10. Differences among the various final state effects are quite small in peripheral collisions as shown in Fig. 10 (a). We see slight depletion in the low p T and enhancement in the high p T regime of the p/π ratio when parton evolutions are included. The impacts of final state rescatterings are more obvious in central collisions shown in Fig. 10 (b). Without any final rescattering effects, the p/π ratio at low p T qualitatively agrees with the data but undershoots in the higher p T region, missing the enhancement feature from peripheral to central collisions suggested by data. The scenario of pure hadron rescatterings as represented by the dotted line shows a large enhancement of the p/π ratio in central pp collisions. The parton rescattering effect seems to be quite similar to that in peripheral collisions. An overall agreement can be reached after both the parton and hadron interactions are included. It is seen in this comparison that both the final state parton and hadron effects are necessary to reproduce the multiplicity dependent baryon to meson enhancement in our transport model approach.

Long range correlation and initial state conditions
Other than the radial flow feature discussed above, particle correlation measurement is a more differential observable to investigate the initial state geometry related collectivity effects in pp collisions. We find that the appearance of the near side ridge structure in two particle long range correlations can be largely related to the initial sub-nucleon fluctuations in transverse space.
In Fig. 11, we present the multiplicity dependence of the projected correlation functions C(∆φ) = 1 Ntrig dN pair d∆φ obtained with the overlap function T (x, y, b) weighting method ( Fig. 11(a)) and the constituent quark assumption ( Fig. 11(b)). The final state parton scattering cross section is 0.2 mb and hadronic rescatterings are switched on in this comparison. The trigger and associate hadrons are selected with the transverse momentum requirement 1 < p T < 3 GeV/c in the acceptance of |η| < 2.4 following the analysis performed at the CMS experiment [8]. Events are separated into two categories based on N sel , the number of selected charge tracks with minimum transverse momentum p T > 0.4 GeV/c within |η| < 2.4. High multiplicity events are required to have N sel > 80 and low multiplicity events are defined with N sel < 20. The two hadrons in each pair must be separated with a pseudo-rapidity gap |∆η| > 2. The correlation function has been corrected following the standard zero-yield-at-minimum procedure. It is suggested by the experimental data that a significant near side ridge structure exists in the correlation function as a local maximum at ∆φ ∼ 0 in the high multiplicity pp events [6,8]. Our study shows a long-range ridge-like structure is present only if the proton matter distribution is modeled based on the constituent quark picture, indicated by the near side peak of black dots in Fig. 11(b). Considering that all the final state interaction parameters are set to be the same in this comparison, this is a quite striking result implying the connection of the induced long range correlation with the underlying sub-nucleon fluctuation effects.
We also ran the same analysis by considering only hadronic rescatterings (0mb parton rescattering cross section) and only parton rescatterings (0.2mb cross section without hFSI) based on the quark constituent assumption in Fig. 11(c) and Fig. 11(d), respectively. The near side ridge persists at high multiplicity when the parton rescattering is included. No near side peak can be found in high multiplicity events with only final hadronic interactions. Therefore, the near side ridge like structure in correlation function is more likely to be developed in the parton evolution stage.
It is interesting to see that the long range near side correlations also exist within the full AMPT model but from a different perspective. Large spatial eccentricity can be induced via the three fire ball like arrangement along the impact parameter direction in AMPT as long as the impact parameter b is not too small to eliminate the separation of different particle sources as already shown in Fig. 4. The spatial profile of the initial partons after string melting from the full AMPT model is presented in Fig. 12. We show the velocity and position of the partons with t p < 5 fm/c at their formation time with impact parameter b=0.1 fm and b=0.6 fm. For central collisions, the partons are roughly placed around the same position. The parton spatial distribution is elongated on the impact parameter direction for peripheral events. Figure 11(b) also shows the long range correlation function from the full AMPT model in the dash lines, wherein the near side ridge arises in high multiplicity pp events due to the parton evolution effects developed from highly eccentric initial geometries [77,50].  We believe that the sizable initial eccentricities created in our current work and the full AMPT model are the key factors to reproduce the flow like long range correlations found in high multiplicity pp data. The separated particle sources in the transverse plane are necessary to generate the large eccentricities for the evolving parton system within the transport model framework. In our current work, the spatial separation is produced through the hot spots like regions due to the sub-nucleon fluctuations using quark constituent assumptions. The full AMPT does not include any sub-nucleon structures but models the initial parton spatial distribution with two or three fire balls along the impact parameter. As the string shoving mechanism also generates significant long range two particle correlations, more detailed studies to understand the implications of these different model implementations are planned in our future work.

Summary
Experimental results revealing collectivity like behaviors in high multiplicity pp collisions have been considered as evidence supporting the creation of deconfined quark gluon medium in small systems. In this work, we introduce a novel transport model approach to systematically study the collective phenomena observed in pp events. We combine PYTHIA8 initial states and AMPT final state interactions together with several options on the proton geometry assumptions in this approach. We show in the study that both parton and hadron final state interactions are important to understand the multiplicity dependent mass ordering in the pp transverse momentum spectra. The near side ridge structure observed in two hadron long range correlations is found to be developed during the partonic phase in our transport model approach. This observation can be regarded as an indication for the creation of deconfined quark matter in high multiplicity pp events. Our study also shows that the appearance of the near side ridge is affected by the proton sub-nucleon fluctuations.
We also note that final state rescattering for pp collisions focused on the hadronic interaction channels has been implemented within the PYTHIA event generator itself [78]. Different space-time structure has been established from the string fragmentation perspective [79], whereas sizable hadronic interaction effects are also observed. It is thus worthwhile to extend the application of our current approach to pA and AA within the Angantyr formalism implemented in PYTHIA8 framework. Examples like this including using UrQMD or the PYTHIA internal hadronic interaction model to handle hadron rescatterings based on the PYTHIA Angantyr space time picture have been extensively applied to describe pA and AA collisions [80,81]. In contrast, our study has the capability of considering parton evolution effects in a coherent picture. The development of this approach will lay a solid foundation for future studies of various mechanisms, such as string shoving and parton/hadron evolutions, within the same model and help us understand the origin of the collective phenomena in pp collisions.