Hadronization of correlated gluon fields

Following an explicit example, we present the chain of steps required for an event-by-event description of hadron production in high energy hadronic and nuclear collisions. We start from incoming nuclei, described in the Color Glass Condensate effective theory, whose collision creates the gluon fields of the glasma. Individual gluons are then sampled from the gluon fields' Husimi (smeared Wigner) distributions, and clustered using a new spacetime based algorithm. Clusters are fed into the Herwig event generator, which performs the hadronization, conserving energy and momentum. We discuss the physical implications of smearing and problems with the quasi particle picture for the studied processes. We compute spectra of charged hadrons and identified particles and their azimuthal momentum anisotropies, and address systematic uncertainties on observables, resulting from the general lack of detailed knowledge of the hadronization mechanism.

Following an explicit example, we present the chain of steps required for an event-by-event description of hadron production in high energy hadronic and nuclear collisions. We start from incoming nuclei, described in the Color Glass Condensate effective theory, whose collision creates the gluon fields of the glasma. Individual gluons are then sampled from the gluon fields' Husimi (smeared Wigner) distributions, and clustered using a new spacetime based algorithm. Clusters are fed into the Herwig event generator, which performs the hadronization, conserving energy and momentum. We discuss the physical implications of smearing and problems with the quasi particle picture for the studied processes. We compute spectra of charged hadrons and identified particles and their azimuthal momentum anisotropies, and address systematic uncertainties on observables, resulting from the general lack of detailed knowledge of the hadronization mechanism.
While the former group employs hydrodynamic calculations, where hadronization is encoded in the equation of state and thus relies on the assumption of thermal equilibrium (with only small deviations due to viscous effects), the latter group has often compared experimental data with parton level results, or employed independent fragmentation, which should only be valid at high p T > 1 − 2 GeV (see e.g. [46] for details). Apart from a study of p+p collisions [40], there has not been an eventby-event calculation in the CGC that employs a sophisticated hadronization description, and would thus allow for direct comparison to experimental data.
In this work, we introduce a new clustering algorithm, dubbed SAHARA (Spacetime Arranging HAdRonizer Application), which can take combined momentum and coordinate space distributions of partons, and generate input for existing hadronization routines, as those available in e.g. PYTHIA [47] or Herwig [48,49]. While the presented framework and clustering algorithm can be applied to any type of partonic description, and in principle be connected to a large variety of hadronization schemes, we will present a specific example consisting of an IP-Glasma initial state calculation in 5.02 TeV p+Pb collisions, determination of the gluon Wigner (Husimi) distribution, clustering with SAHARA, and hadroniza-tion of individual clusters using Herwig.
We compute charged hadron transverse momentum spectra and azimuthal anisotropies v n and v n (p T ), as well as v n (p T ) for identified particles. We further study the dependence of calculated observables on the clustering parameters and compare parton to hadron level results. Our results provide a first direct comparison of event-byevent CGC calculations with experimental data in p+Pb collisions, when no intermediate hydrodynamic stage is included.
This paper is organized as follows: In Section II we describe how initial gluon distributions are obtained from solving the Yang-Mills equations, computing the gluon Wigner distributions, and smearing and sampling the distributions. Section III introduces the clustering procedure that results in clusters of gluons that will be hadronized independently. The hadronization is done via Herwig, which is briefly described in Section IV. We present results for observables in Section V and discuss the dependence on clustering parameters. We conclude in Section VI. We finally present distributions of the invariant masses of SAHARA and Herwig clusters in Appendix A.

II. GLUON HUSIMI DISTRIBUTION FROM THE IP-GLASMA
The IP-Glasma model is an event-by-event implementation of the leading order CGC framework that computes the gluon fields produced in a high energy nuclear collision and their time evolution based on solutions of the classical Yang-Mills equations [50,51]. The energy momentum tensor of these gluon fields has been used as input for event-by-event hydrodynamic calculations whose results have shown very good agreement with experimental data for a large variety of collision systems (see e.g. the recent comprehensive study in [52]). Here, we extract the (boost invariant) gluon Wigner distribution as done previously in [53]. This involves evaluating equal time correlation functions in Coulomb gauge and the projection onto transverse polarization states of the free theory [54] dN Wigner which we compute at the time τ = 0.2 fm/c, after which the system becomes essentially free-streaming [38,51]. Here, λ runs over the transverse polarizations, a over the N 2 c − 1 colors. The time dependent transverse polarization vectors are given by ξ (λ) p T ,µ (τ ), and in Coulomb gauge take the form [38,54] with p T = (p x , p y ) and H (2) α the Hankel functions of the second kind and order α.
The gluon fields as a function of transverse position are given by A a µ (x T ) and obtained from the IP-Glasma calculation, described e.g. in [38,[50][51][52]. The model includes fluctuations of nucleon positions and three subnucluonic hot-spots, whose distributions, along with the IP-Sat model [55,56] that provides the saturation scale Q s for a given collision energy, rapidity, and thickness function, determine the spatial color charge density distributions in the incoming proton and nucleus. The setup is as described in [52], including normalization fluctuations in the Q s of each hot spot, except that some parameters are chosen differently. In particular, here we choose the infrared regulator m = 0.4 GeV, and the width parameter of the normalization fluctuation σ = 0.5. Using the resulting color charge density, the McLerran-Venugopalan model [57,58] determines the fluctuating color charges, which form the external current in the Yang-Mills equations for the incoming gluon fields [59]. We then numerically solve for the gluon fields produced in the collision and evolve them with the source free Yang-Mills equations to τ = 0.2 fm, which corresponds to approximately a time scale of 1/Q s .
We stress that although the Wigner distribution (1) contains all information about single particle states, it is not positive semi-definite in all phase-space regions. Consequently, it is not a probability distribution, which is necessary for a quasi-particle interpretation. Since in proton-proton or proton-nucleus collisions the spatial variations on scales of the order of the proton size R p ∼ 0.4 fm, occur on essentially the same scale as the de-Broglie wave-length of typical excitations λ Qs ∼ 0.2 fm, a quasi-particle interpretation is problematic, as there is no clear separation of scales between the length scale of gradients and the quantum mechanical size of the wavepacket of a single particle. Hence, in order to obtain a positive definite quasi-particle distribution, it is necessary to perform a coarse graining procedure before the sampling of individual gluons, which then can be clustered and hadronized.
In order to obtain a quasi-particle distribution from the Wigner distribution (1), we need to smear it over phasespace volumes of size σ x σ p ≥ . This leads to a single particle distribution, the so called Husimi distribution [60], which within our boost-invariant picture reads We use σ x = 0.197 fm and σ p = 1 GeV, to achieve a reasonable compromise between spatial and momentum resolution. Since the classical Yang-Mills calculation yields results for g 2 dNg dηsd 2x T d 2p T rather than the multiplicity dNg dηsd 2x T d 2p T , the strong coupling constant g can be adjusted a posteriori to produce the correct charged particle multiplicities after hadronization, resulting in g = 2.75 for our study. The various steps in this procedure are illustrated in Fig. 1, where in the different panels we present the spatial distribution of the energy density per unit rapidity g 2 τ (x, y), the momentum space distributions g 2 dN g /dydp x dp y , as well as the Wigner and Husimi distributions. We note that for the energy density, the result before smearing is obtained from the gauge invariant operator definition, while after the smearing the local energy density is calculated from the momentum integral of the phase-space distribution in Eq. (4). Generally, one observes that the minimal uncertainty smearing has a non-negligible effect on both the spatial and momentum distributions, indicating that any theoretical description based on localized quasi-particles is definitely pushed to its limits of validity in p+Pb collisions. Nevertheless, even though fluctuations on short distance and momentum scales are completely washed out by the smearing, the structure of the event on distance scales σ x and momentum scales σ p remains intact.

III. CLUSTERING WITH SAHARA
Since the hadronization in typical high energy physics event generators requires the notion of events with indi-  Figure 1. Left: Spatial profile of the energy density per unit rapidity g 2 τ (x, y) in the transverse plane, before (top) and after (bottom) the minimal uncertainty smearing. Center: Momentum space distribution of gluons g 2 dNg/dyd 2 pT in the transverse plane, before (top) and after (bottom) the minimal uncertainty smearing. Right: Wigner distribution (top) and Husimi distribution (bottom) as a function of transverse momentum, evaluated at 0.08 fm < x < 0.14 fm and −0.24 fm < y < −0.18 fm.
vidual partons, we subsequently sample a collection of individual gluons over a spacetime rapidity range determined by η max = 2, with the transverse coordinates x T and p T assigned according to the Husimi distribution for a given event in Eq. (4). Evidently, a single sampling of the event will not capture all the features of the underlying IP-Glasma event. Therefore, we sample each IP-Glasma event multiple times, and perform the subsequent hadronization procedure independently for each sample.
Starting from an individual event, characterized by the spacetime positions (t g i , x g i ) and four momenta (E g i , p g i ) (with i = 1 . . . N g ) of the produced gluons, we invoke SAHARA to arrange the event into color neutral clusters, which will hadronize independently of each other. Since the gluons initially produced in a (semi-) hard scattering can have undergone additional re-interactions in the final state, the clustering in SAHARA is not tied to the hard production process, but instead based on the concept of spacetime locality of the hadronization process, where gluons i, j with a small distance of closest approach in the center of mass frame of the collision are likely to hadronize together, as they have closely encountered each other over the course of the spacetime evolution of the collision. We further note that by defining the distance in Eq. (5) in terms of the (forward restricted t > 0) distance of closest approach, any free-streaming evolution in the final state does not change the proximity measure and the cluster hadronization will be modified only if the produced gluons experience final state interactions. Next, in order to implement the idea of a local hadronization process, we define a cluster action S (c) cl and perform a statistical clustering, where the probability to obtain a given cluster configuration {C} with N cl nonempty clusters is given by with Z = {C} P ({C}) being the partition function. Specifically, for each cluster c containing i, j = 1, · · · , N g c gluons, we define Cluster 0 Cluster 1 Cluster 2 Cluster 3 Cluster 4 Cluster 5 Cluster 0 Cluster 1 Cluster 2 Cluster 3 Cluster 4 Cluster 5 Cluster 0 Cluster 1 Cluster 2 Cluster 3 Cluster 4 Cluster 5 Cluster 0 Cluster 1 Cluster 2 Cluster 3 Cluster 4 Cluster 5 Figure 2. Example of a cluster configuration in SAHARA for a typical p+Pb event with Ng ≈ Ng , with α = 0 (top) and α = 10 (bottom). Shown is the three-dimensional x, y, η structure (left) and the projection onto the transverse x,y, plane (right). Different SAHARA clusters are indicated by different colors and contours on the right show the smeared energy density distribution.
] is the average distance of closest approach of gluons inside the cluster and N g c is the number of gluons inside the cluster, such that the factor N g c d2 c can be thought of as the average length of the color flux tube connecting individual gluons. On the other hand, the clustering parameterᾱ = α fm −1 controls the importance of spacetime locality within each cluster, relative to the combinatorial possibilities of clustering. Evidently, to achieve color neutral clusters, each cluster needs to contain at least N min = 2 gluons. However, instead of requiring a minimal number of two gluons per cluster, we will treat N min ≥ 2 as a free parameter and vary its value to study the effect on hadronic observables.
Statistical clustering in SAHARA is performed by a Markov Chain Monte Carlo procedure, which consists of three distinct steps, where clusters can split, merge or exchange particles to ensure ergodicity. Since the statistical clustering is based on the probabilities in Eq. (6) it is independent of the details of the algorithmic implementation and we refrain from providing an exhaustive discussion of technical details.
We present an example of a SAHARA cluster configu-ration in Fig. 2, for the same 5.02 TeV p+Pb event with N g ≈ N g ≡ N minbias (where · is the average over all events) depicted in Fig. 1. Points mark the position of gluons and arrows indicate the corresponding velocities.
Different SAHARA clusters (0,. . . ,5) are indicated by different color coding. On the left, we show the three dimensional distribution of the gluons' spatial positions. Upper panels correspond to a purely statistical clustering (α = 0), and lower panels to clustering with preferably small distance of closest approach (α = 10). In that case, one can see that clusters extend over several units in rapidity, but are well localized in the transverse x-y plane. On the right hand side we show projections to the x-y plane overlayed with the smeared energy density distribution of the event at τ = 0.2 fm/c. The stronger localization of clusters for the α = 10 case (bottom) compared to α = 0 (top) is clearly visible here.

IV. HADRONIZATION WITH HERWIG
Each individual cluster generated by SAHARA is then hadronized using Herwig 7 [48,49]. In order to do so, we need to assign color connections to the gluons in each cluster. Since we have sampled individual gluons from the single particle distribution (4), we have no information on the color flow of the event or individual clusters, and we therefore choose to perform a democratic color assignment in the spirit of maximizing the associated entropy. In practice, we need to assign two "color values" to every gluon, which indicates its connection to one or two other gluons (we do not include quarks or anti-quarks, which would have a single color index). Note that one in principle could employ a different prescription, e.g. connecting gluons according to their invariant mass or another measure inspired by the kinematic dependence of color subamplitudes. Since the masses of the SAHARA clusters are already peaked at small values (see Appendix A) and do not constitute a genuine hard jet configuration, we do not consider any more sophisticated treatment in this context. We have also found that color reconnection in Herwig (see [61] and references therein) within the SAHARA clusters do not affect our results and that we obtain a reasonable mass spectrum of Herwig clusters, which supports our choice. Because the total event record, i.e., the cluster to be hadronized, has to be color neutral, we begin with "color neutral gluons", i.e., gluons that have the same color and anti-color index. We then proceed to independently shuffle the color and anti-color values of all gluons until no color neutral gluons are left in the final event record.
For every cluster we provide Herwig with a Les Houches [62] event record that contains the parton level information, including the gluons' color connections and four-momenta. While Herwig can in principle use coordinate-space information for its color reconnection, we here do not include this, because the SAHARA clustering already provides the necessary correlations in color space. In practice, Herwig is provided with a singlet-intogluons type of process in the center-of-mass frame of each individual SAHARA cluster.
As some of the SAHARA clusters can still have significant masses, and no additional parton cascade is included in the simulation, we do let Herwig perform the parton showering according to the coherent branching algorithm outlined in [63]. After the gluons have been split into quark-anti-quark pairs, clusters of color connected quarks and anti-quarks are formed. 1 Their mass spectrum is compatible with the assumptions of the cluster hadronization model [64], and mostly independent of the choice of α (see Appendix A), and they undergo the usual mechanism of cluster fission and cluster decay into hadrons. Given the universality of the cluster hadronization model in dependence on the center-of-mass energy, provided that there has been a coherent branching evolution, our approach should be viable for the hadronization 1 We stress that one should not confuse the clusters generated by SAHARA with the ones that Herwig itself produces within a SAHARA cluster. We discuss the invariant mass distributions of SAHARA clusters and Herwig clusters in Appendix A.
of gluon clusters originating from a glasma simulation.

V. RESULTS
The Herwig output consists of final particles and their momenta after showers and resonance decays. Combining results from all clusters and applying the appropriate boosts from the cluster's center of momentum to the lab frame, we obtain an event-by-event record of the produced hadrons. We oversample a single IP-Glasma event several thousand times, meaning that we sample individual gluons and run SAHARA followed by Herwig many times, and combine all oversampled results to obtain the final particle spectra for each IP-Glasma event.
We analyze those using the publicly available toolkit by Chun Shen [65], which has been used for a variety of hybrid hydrodynamic and hadronic cascade calculations [16,52,66,67], and obtain transverse momentum spectra and azimuthal momentum anisotropies in different centrality classes. We will compare hadron level results to those for gluons and study the dependence on the clustering parameters α and N min .
We begin our analysis by comparing the final charged hadron transverse momentum spectrum to the gluon spectra before and after smearing in the centrality class around twice the minimum bias multiplicity N minbias in √ s = 5.02 TeV p+Pb collisions in Fig. 3. We also compare to the hadron spectrum obtained from performing independent fragmentation of the smeared gluon spectrum using the next to leading order (NLO) Kniehl-Kramer-Potter (KKP) fragmentation functions [46] as described in [68].
We find that for the largest studied transverse momenta, 5 GeV < p T < 8 GeV, the result from IP-Glasma+SAHARA+Herwig agrees well with the IP-Glasma+KKP result, which can be expected. For decreasing transverse momenta, differences increase, with the largest discrepancy for p T 0.5 GeV, where the application of fragmentation functions is indeed questionable.
In comparison to the gluon spectra we observe that the hadron spectrum is significantly steeper, as expected from earlier studies [68]. We also note that the effect of smearing on the gluon spectrum is rather mild, only becoming significant for p T 1 GeV.
We finally show the experimental data for the non single diffractive (NSD) charged hadron spectrum measured by the ALICE Collaboration [69] and scaled by a factor of 2, as we show our 1.95 < N < 2.05 bin. This is not a perfect comparison, as the ALICE result contains all NSD events and not just a small selection like our result, which may affect the shape of the p T spectrum. Nevertheless, the comparison shows rather good agreement, particularly at the lowest and highest p T shown, with small deviations in the range 2 GeV < p T < 6 GeV.
We note that in Fig. 3, we used the clustering parameter α = 10 and minimal number of gluons per cluster Hadron spectra from the IP-Glasma+SAHARA+Herwig calculation (thick solid line) compared to hadrons obtained from hadronizing IP-Glasma gluons using NLO KKP fragmentation functions (thin solid line). We also show gluon spectra from IP-Glasma (dotted line) and the gluon distribution after smearing (dashed line). We compare to the non single diffractive charged hadron spectrum measured by the ALICE Collaboration [69], scaled by a factor of 2, as we show the bin around 2N minbias .
N min = 4. Changing α to 0 resulted in a 5-10% difference in the charged hadron p T -spectrum and we will return to the dependence on clustering parameters below, where we discuss azimuthal anisotropies, which are more sensitive to α.
Next, we will analyze the azimuthal anisotropies of produced particles, which are determined using the scalar product method [70,71] where the anisotropic flow vectors are given by and The angles φ j are the azimuthal angles of the j'th particle. The sum over j runs over all final particles in either a specific p T bin and 0.5 < η < 2, after oversampling, where η is pseudo-rapidity, or in the case of Q ref n in a reference bin, chosen to cover −2 < η < −0.5 and 0.2 GeV < p T < 3 GeV. The number of particles in the reference bin is N ref and which has self correlations removed by means of the second term in the numerator.  Note that in the case of gluons, we use the same ranges for the reference (and p T -) bins. For the analysis of identified particles, we constrain particles in the p T -bin to be of the desired species, and correlate them with all charged hadrons in the reference bin. For integrated v n {2} = C n {2} 1/2 , where C n {2} is defined in analogy to C ref n {2}, but using different rapidity intervals for Q * n and Q n , respectively, introducing a gap to eliminate residual non-flow (which is already reduced by oversampling). The subtraction of self correlations is unnecessary in that case.
Results for v 2 {2}(p T ) of charged hadrons and gluons are shown in Fig. 4 for events with approximately two times the minimum bias multiplicity N = 1.95 − 2.05 N . Solid lines in Fig. 4 correspond to the final result for charged hadrons, which we compare to the gluon v 2 {2}(p T ) shown by filled circles. It is interesting to note that the maximal value of the charged hadron v 2 {2}(p T ) is comparable to that of the gluons. Furthermore, the charged hadron v 2 {2}(p T ) is larger than that of the gluons at low p T , which is expected, as during fragmentation larger p T gluons fragment into smaller p T hadrons. The larger hadron v 2 {2}(p T ) compared to gluons for p T > 4 GeV could be explained by realizing that azimuthal anisotropies of gluons from the IP-Glasma decorrelate quickly when their p T difference is increased [38], while the additional "smearing" in p T from hadronization reduces this effect. Thus, when p T moves out of the range of the reference bin, 0.2 GeV < p T < 3 GeV, the decorrelation is more rapid for gluons than charged hadrons.
We note that the charged hadron v 2 {2}(p T ) in the experimental data is approximately a factor of three larger (for e.g. p T = 2 GeV) than our result [72,73]. The discrepancy likely stems from the presence of large final state effects in nature that are neglected here and can e.g. be implemented using a partonic transport approach [53] or by incorporating a hydrodynamic phase [52].
We also show the comparison to the charged hadron v 2 {2}(p T ) that one obtains by simply folding the gluon distribution with the KKP fragmentation functions. In this case, the p T dependence is very different from our result, as the fragmentation functions mostly lead to a shift of the maximal v 2 to lower momenta, along with a smearing in p T that reduces the overall magnitude. We note that in [32] different fragmentation functions were explored, and a strong dependence on the near side associated yields (equivalently, the v 2 ) was found. Since the azimuthal anistotropies of charged particles are sensitive to non-perturbative hadronization effects, a realistic modeling is required to compare theoretical results with experimental data. Within our SAHARA+Herwig framework, we have also investigated the effect of disabling color reconnections in Herwig and found no effect on the results for v 2 (p T ).
In order to understand if an artificial anisotropy can be generated in the hadronization process, and to estimate its size in comparison to the full result, we ran the same simulation after randomizing the azimuthal angles of all gluons in an event. The result is shown in Fig. 5, together with the full result, that we repeat from Fig. 4. The expected zero v 2 {2}(p T ) for gluons is verified. For hadrons, however, a finite v 2 {2}(p T ) is observed (dashed line in Fig. 5), albeit significantly smaller than that in the full simulation that includes the gluon momentum anisotropy. The likely origin of this anisotropy for hadrons are the changing off-diagonal elements of the energy momentum tensor 2 during the hadronization pro-2 While energy and momentum are conserved, e.g. T xy is not cess, which, as discussed in [74], can induce azimuthal anisotropies in momentum space.
The third harmonic v 3 {2}(p T ) for charged hadrons and gluons in the 2N minbias multiplicity class is shown in Fig. 6. We find that the charged hadron v 3 {2}(p T ) resembles that for gluons closely for all p T . Again, experimental data for v 3 {2}(p T ) is approximately three times larger than our result [72,73], as we lack strong final state interactions. Again, we show for comparison the result obtained when folding the gluon distribution with NLO KKP fragmentation functions. As in the case of v 2 {2}(p T ), the transverse momentum shape of v 3 {2}(p T ) is very different in this case, with the largest values reached at much lower p T .
As discussed previously in the literature [17,38], the momentum anisotropy of gluons in the CGC decreases as a function of multiplicity, which is intuitively explained within the color domain model [36,75]. In Fig. 7 we show the charged hadron v 2 {2}(p T ) for three different multiplicity classes. While statistical errors are large, we see a systematic decrease of v 2 {2}(p T ) in the bin containing events with multiplicity > 3.45 N minbias compared to the 2N minbias class. The difference between the classes around N minbias and 2N minbias is small. Fig. 8 shows the v 2 {2}(p T ) for identified pions, kaons, and (anti-)protons. While statistical errors are large, one can see a clear mass ordering of the v 2 {2}(p T ). This is expected, as particles emerge from clusters with a common transverse velocity, similar to the case of particlization from a fluid cell in a hydrodynamic simulation. The effect of mass splitting was also observed in calculations using IP-Glasma and PYTHIA for p+p collisions for similar reasons [40]. Given these observations, it is hard to conserved during our procedure.  imagine a hadronization mechanism that does not lead to mass splitting of the v n , and identifying the correct underlying mechanism requires a detailed quantitative analysis.
So far we have analyzed azimuthal correlations for the default SAHARA clustering with α = 10 and N min = 4. Next, to estimate the uncertainty from our clustering algorithm, we will vary the parameters α and N min discussed below Eq. (7). Fig. 9 shows v 2 {2}(p T ) in the 2N minbias multiplicity class for our standard choice of α = 10 together with the case α = 0. As vanishing α means that gluons are randomly assigned to clusters, independently from their position and momenta, it is not surprising that the resulting v 2 {2}(p T ) is reduced in this case. Increasing α prefers clustering of nearby gluons (defined via the distance of closest approach), which leads to larger cluster momenta, which in turn will better retain the azimuthal anisotropy of the gluons.
In Fig. 10 we study the dependence of charged hadron v 2 {2}(p T ) on the minimal number of gluons per cluster N min . The results for N min = 3, 4 and 6 agree within our statistical errors, but we see that the case N min = 4 leads to systematically larger v 2 {2}(p T ) than the other two. While it is expected that an increasing N min leads to a reduction of v 2 {2}(p T ), because clusters will be forced to contain more gluons with increasingly different positions and momenta, it is not obvious why reducing N min from four to three reduces v 2 {2}(p T ). One possible explanation is that clusters with fewer particles also have a smaller extent in rapidity, such that correlations within individual clusters are less likely to survive the rapidity gap introduced in Eq. (8). However, other effects associated e.g. to details of the color structure, caused by having fewer particles in a cluster, are also conceivable.
We continue with the study of charged hadron mean  Figure 11. The charged hadron mean transverse momentum pT as a function of charged hadron multiplicity (computed in three different bins) compared to experimental data from the ALICE Collaboration [76]. For the 2N minbias bin we show results for different choices of α and Nmin.
transverse momentum p T as a function of charged hadron multiplicity in Fig. 11. We compare to experimental data from the ALICE Collaboration [76], and for the 2N minbias centrality class vary α and N min to estimate the uncertainty from the clustering procedure. Generally, we find that the agreement with the experimental data is rather good, in particular for the lower multiplicities, while the dependence on the clustering parameters α and N min in the studied range is weak.
In Fig. 12 we show the p T -integrated charged hadron v n {2} for n = 2, 3 as a function of charged hadron multiplicity. We do not show a comparison to experimental data, as it is at least a factor of three above our results, as discussed above. Similar to Fig. 11, we again show results obtained for varying α and N min for the 2N minbias bin. We find that for the integrated v 2 {2} the dependence on the parameters is rather small as it is dominated primarily by low p T particles, for which we had seen a small effect from the choice of α or N min . Even though for v 3 {2} the dependence on the parameters is slightly larger, it is hard to deduce a clear systematic effect. Besides the results from full IP-Glasma+SAHARA+Herwig simulations, which contain initial state momentum correlations of gluons in IP-Glasma, we also show results for the case of randomized gluon azimuthal angles. While the charged hadron v 3 {2} is consistent with zero in this case, v 2 {2} is finite, as we had observed for v 2 {2}(p T ).
Finally, we present the ratios of kaons and protons to pions as functions of charged hadron multiplicities in Fig. 13. Somewhat surprisingly, as we do not include any quark degrees of freedom, the ratio of kaons to pions agrees well with the experimental data from the ALICE Collaboration [77]. The ratio of protons to pions underestimates the experimental data by approximately 30-40%. We also show results for α = 0, which are very similar to those for α = 10, indicating that the hadrochemistry is  Figure 13. Particle ratios K/π and p/π as functions of charged hadron multiplicity compared to experimental data from the ALICE Collaboration [77]. not significantly affected by the details of the clustering.

VI. CONCLUSIONS
We presented a newly developed hadronization scheme based on spacetime locality along with a first application to gluonic systems produced in p+Pb collisions at the LHC. As the gluon Wigner distribution is not positive semi-definite, a quasi-particle interpretation is not immediately possible. This is in principle a fundamental problem, which we circumvent by smearing the Wigner distribution over phase-space volumes determined by the uncertainty principle, which makes it positive semi-definite and allows the sampling of individual gluons. These gluons are then clustered based on their distance of closest approach, and the clusters are passed to Herwig, which performs the hadronization according to its default procedure.
We computed transverse momentum spectra and azimuthal momentum anisotropies of charged hadrons, showing that the hadronic observables v n (p T ) closely resemble the gluonic ones, with the largest differences for v 2 (p T ) at the lowest p T 1.5 GeV. This validates to some degree the comparisons of parton level results to experimental data, as done with some previous CGC based calculations (e.g. [43,44]). Performing independent fragmentation by folding gluon distributions with NLO KKP fragmentation functions on the other hand modifies the p T dependence of the v n (p T ) dramatically.
While spectra and mean transverse momenta agree well with the experimental data, the overall magnitude of the v n is significantly smaller than what has been observed in experiment, which supports the conclusion [16,78,79] that final state effects (which we do not include here) are crucial to correctly describe the azimuthal anisotropies, also in small systems.
For identified particles we observe a mass splitting of v 2 (p T ). It is qualitatively similar to what hydrodynamic calculations predict and what was observed in the experimental data. Further, the particle ratios of K/π and p/π are reasonably well described, given that we do not include any quark degrees of freedom.
We also studied the dependence of the clustering parameters on the final observables, finding only a relatively weak dependence, highlighting the robustness of the prescription. Since the SAHARA framework provides a general interface to the hadronization mechanisms implemented in high-energy physics event generators, it would be interesting in the future to also couple to other hadronization descriptions, such as that used in PYTHIA.
Evidently, to achieve a satisfactory description of collective phenomena in small systems, it will be important to include final state effects within our framework. One possible direction would be to follow [53] and introduce a phase of partonic scatterings via solutions of the Boltzmann equation (for example using BAMPS [80]), whose outcome can then be clustered using SAHARA and hadronized. This would provide a complete microscopic description of nuclear collisions from the beginning to end and could be applied to a plethora of processes.

VII. ACKNOWLEDGMENTS.
We thank Chun Shen for useful discussions and advice on using the iEbE* framework. We thank Prithwish Tribedy and Raju Venugopalan for helpful discussions. This work was supported by the Helmholtz International Center for FAIR within the framework of the LOEWE program launched by the State of Hesse. This work is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the CRC-TR 211 'Strong-interaction matter under ex-