Hadronization Scheme Dependence of Long-Range Azimuthal Harmonics in High Energy p+A Reactions

We compare the distortion effects of three popular final-state hadronization schemes. We show how hadronization modifies the initial-state gluon correlations in high energy p+A collisions. The three models considered are (1) LPH: local parton-hadron duality, (2) CPR: collinear parton- hadron resonance independent fragmentation, and (3) LUND: color string hadronization. The strong initial-state azimuthal asymmetries are generated using the GLVB model for non-abelian gluon bremsstrahlung, assuming a saturation scale Qsat = 2 GeV. Long-range elliptic and triangular harmonics for the final hadron pairs are compared based on the three hadronization schemes. Our analysis shows that the process of hadronization causes major distortions of the partonic azimuthal harmonics for transverse momenta at least up to pT = 3GeV. In particular, they appear to be greatly reduced for pT<1{\div}2GeV.


Introduction and motivation
Multi-particle, long-range in pseudo-rapidity azimuthal correlations are widely studied in high-energy nuclear collisions at RHIC Au+Au [1,2] and LHC Pb+Pb [3][4][5]. In particular, they are considered as a signature for the "perfect fluid" behavior of the strongly coupled Quark-Gluon-Plasma (sQGP) produced in such reactions.
The recent discovery of long-range p+A azimuthal harmonics , see Fig. 1(a), with magnitudes and p T dependence comparable to the ones found in A+A collisions [6][7][8][9], and the near energy independence of these A+A moments observed in the Beam Energy Scan (BES) at RHIC [10] together with LHC have challenged the uniqueness of the sQGP interpretation of v n in A+A. Prior to the recent p+A and BES data, it was assumed that in smaller transverse size p+A systems or lower energy A+A reactions, the perfect fluid sQGP "core" would gradually change into a highly dissipative hadron resonance gas "corona" and lead to a substantially smaller magnitude of the azimuthal harmonic moments. The observed near independence of v n (p T ) to system size and initial energy density has motivated the search for possible alternative non-hydrodynamic sources of azimuthal harmonics. Experimental v n harmonics for p+Pb collisions and |∆η| > 2 taken from ATLAS [12]. The magnitudes of the moments are comparable to the one observed in A+A reactions. Right panel: Azimuthal Fourier moments arising from the GB distribution (Eq. (1)) as a function of the gluon transverse momentum and for different choices of the exchanged momentum, as computed in [22].
While it has been shown that hydrodynamic equations, with particular assumptions about the initial and freeze-out conditions [11], are sufficient to describe the data, the uniqueness of that description is not obvious, especially given the unexpected features of the previously mentioned data. One important class of non-hydrodynamic models proposed to explain this puzzle is based on the Color Glass Condensate (CGC) and Glasma paradigm involving initial-state nonperturbative classical field correlations controlled by a gluon saturation scale, Q sat (x, A) [13][14][15][16][17][18][19][20][21]. A simpler perturbative QCD source of multi-gluon azimuthal correlations due to classical nonabelian field interference effect was recently presented in GLVB [22] based on the well-known Gunion-Bertsch (GB) LO analytic formula [23]. This non-abelian gluon radiation generates nonzero v n moments with a shape that closely resembles the experimental one -see Fig. 1(b). The magnitudes of those parton level harmonics are however too large and some kind of damping mechanism is therefore required. A natural candidate for this job is hadronization. In fact, the primary advantage of GLVB, for our purpose of exploring systematic theoretical uncertainties associated with hadronization scheme dependence of azimuthal harmonics, is its ease of adaptability to Monte Carlo (MC) multi-particle production in p+p, p+A and A+A collisions via the HIJING [24] type of models as emphasized in [22].
Thus far CGC/Glasma phenomenology of azimuthal harmonics moments in p+A and A+A reactions has been limited to the idealized local parton-hadron hadronization scheme [25] that by assumption preserves the parton level correlations into the final-state hadrons. This guarantees minimal distortion of the initial-state multi-parton correlations predicted by QCD. More realistically, it is well-known from decades of hadronic phenomenology that phenomena like the production and dacay of intermediate hadron resonances and final-state correlations can introduce non-trivial complications and, therefore, some non-perturbative hadronization scheme is required. Hadronization phenomenology cannot be rigorously predicted from QCD, but two generic approaches, carefully tuned to e + e − , ep, and pp data, have been developed during the 2 years. These models are the independent fragmentation scheme [26] and the LUND string model [27] implemented in the JETSET algorithm [28]. They can be used to estimate the distortions of initial-state partonic correlations due to a more realistic hadronization process. In particular, in the nuclear event generator HIJING, JETSET is used to convert the multiple diquarkquark beam jets with gluon mini-jet kinks into the multi-hadron resonances with particle data book properties and decay branchings.
In this paper we study the hadronization scheme dependence of final-state correlations using three different hadronization models that can be conveniently selected within the JETSET code. We deliberately neglect any intermediate parton or hadron level transport or hydrodynamic effect to concentrate exclusively on final-state hadronization modification of initial multi-gluon correlations. We compare: (1) local parton-hadron duality (LPH), (2) collinear (independent) partonto-hadron resonance fragmentation/decay (CPR) and (3) LUND string hadronization models. We treat Q sat as an input parameter controlling the p T range of azimuthal asymmetry of Gunion-Bertsch multi-beam gluon bremsstrahlung. Our results strongly support the analysis of Skokov et al. [21], indicating that the CGC/Glasma predictions (without detailed hadronization modelling) should be limited to the kinematic range p T > 3 GeV to avoid complications and theoretical uncertainty of non-perturbative hadronization physics.
At this point, it should be stressed that the configurations that we will consider in this work are quite simple and hence cannot be considered as a realistic representation of the ridge effect observed experimentally in p+A reactions. However, as we will see, they are already enough to draw strong conclusions about the hadronization mechanism dependence of the final azimuthal harmonics.

Elliptic interference harmonics of gluons from two recoiling beam jets
For our simulations, we used the MC hadronization algorithm JETSET 7.4 [28] with 30 × 10 6 simulated GLVB "events" with two recoiling beam jets. Each beam jet is represented by a high invariant mass qq pair along a "beam axis",ẑ. For our simulations the invariant mass of each beam jet is taken to be 100 GeV. The two beam jets are assumed to scatter with equal but opposite net momentum transfer Q 1 = −Q 2 = (q, ψ), with magnitude distributed according to a simple Rutherford form Q 2 sat /(q 2 + Q 2 sat ) 2 , and azimuthal direction, ψ, distributed uniformly in [−π, π]. For our numerical simulations, we further assumed Poisson fluctuations of the number, N g , of bremsstrahlung gluons with N g = 8 and distributed uniformly in rapidity between kinematic bounds.
In the LUND color string model [27], all gluons are represented by "kinks" of the qq string that deform it in the transverse plane. In contrast, from an independent fragmentation point of view [26], gluons are assumed to be just isolated partons that hadronize independently from each other. We generate the bremsstrahlung gluons transverse momenta and azimuthal angles, k = (k, φ), from the perturbative regularized GB distribution: with Q sat = 2 GeV being the typical momentum scale expected from the CGC model for p+A collisions. We take µ = 300 MeV and Λ QCD = 200 MeV as infrared regulating scales. This setup produces initial azimuthal anisotropy at the parton level that with the simplest LPH scheme leads to harmonic moments similar to experiment, as seen in Fig.1. In particular, 3 φ ∆ the pseudo-rapidity independence and the preference of radiated gluons transverse momenta k to be aligned along ±Q gives rise to a "ridge-like" structure with azimuthal asymmetry resembling the one observed experimentally.
After the hadronization, an analysis is performed over pairs of final pions with criteria as close as possible to [8]. In particular, our pseudo-rapidity cuts on the final pions are |η| < 2.5 and |∆η| > 2 (the so-called long-range), ∆η = η a − η b being the relative psuedo-rapidity of the pair. Both pions are taken in the same p T bin and the correlation functions we computed are given by: where ∆φ = φ a − φ b is the relative azimuthal angle between the two pions and S and B represent the same and mixed event pair distributions respectively [8]. The analysis is repeated for the three hadronization schemes to reveal the hadronization scheme dependence of final azimuthal harmonics shown in Figs. 2 and 3. From these distributions we extracted the n ≤ 6 single-particle moments, v n (p T ), using the following fitting function: The results of this fit for n ≤ 4 are shown in Fig. 4.

Numerical Results
The computation of the two-particle ∆φ distributions as reported in Figs. 2 and 3 shows that the process of hadronization has some deeply non-trivial consequences for the correlations between final state particles. The gluon curves clearly present the behavior expected from the discussion around Eq. (1). The long-range near-side (∆φ 0) peak is simply due to the fact that the gluons are produced with transverse momentum preferentially close to the exchanged momentum Q, while the awayside (∆φ π) peak is produced by the recoil of the two strings. These features extend over the full p T range.
The distributions of the final pairs of hadrons, instead, have some very peculiar properties, depending on which p T region we are considering. For small values of the pion transverse momentum -smaller than a scale, λ, to be determined in Sec. 3 -the initial anisotropy is extremely reduced and the two-pion distributions become more and more uniform for decreasing p T . This behavior is common to both the LUND model and the independent fragmentation model, even though the shape of the ∆φ distributions already appears to be different for p T > 1 GeV. The reasons for this strong dilution of the initial anisotropy are proper to the process of hadronization itself and will be analyzed again in Sec. 3. If a simple system like ours already has such a "decoherence power", one can reasonably expect the same thing to be valid for a more realistic situation -see for example Sec. 4. It should also be stressed that, since the low-p T pions seem to carry no information about the initial gluon distribution, this feature appears be true in general, no matter what causes the partonic correlation in the first place, and hence can be applied to other initial-state models as well.
The second result refers to the higher transverse momentum regime (p T λ). In both models this region presents strong two-particle correlations among the final pions, since the gluon anisotropy is transmitted more efficiently. However, one can immediately appreciate an essential difference between the independent fragmentation and the QCD string models. In the former, by definition, each parton fragments independently from the others and hence the pion correlation function closely resembles the gluon one since, when higher p T are considered, no other relevant sources of angular correlation come into play. For the case of the LUND model, instead, while the near-side peak is generally lower for the pions than for the gluons, an extremely pronounced away-side peak is present. This large ∆φ π signal (several times bigger than the ∆φ 0 one) is a direct consequence of the transverse momentum conservation taking place at each breaking of the color string [28]. Each qq system, in fact, has zero total p T and the transverse fluctuations are governed by a roughly Gaussian probability distribution. This means that when the string breaks and a new pair of partons, say q andq , is created from the vacuum, the p T of these partons has a probability distribution given by with m being the mass of the produced q and κ 1 GeV/fm being the string constant. By conservation of total momentum they are always produced back-to-back. It then follows that, for every hadron of momentum p T there will always be another with momentum −p T , i.e. such that ∆φ = π. This essentially means that the differential distribution C(∆φ; p T ) contains a contribution proportional to δ(∆φ − π), which gets broader when finite p T bins are taken into account, as in any realistic situation. We conclude that the QCD strings hadronization scheme that conserves event energy-momentum (unlike independent fragmentation) can introduce large away-side correlations in the spectrum of the final hadrons that can lead to negative odd harmonics absent at the purely partonic level.  Table 1: Ratios between the pion even moments (v π n ) and the gluon ones (v g n ) for p T < 2 GeV. One immediately notices that there is a consistent damping of the initial magnitudes. In particular, the final-state v n are essentially zero (less than 30% of the initial ones) for p T < 1 GeV. The suppression for p T > 2 GeV appears to be negligible.
All these properties lead to the azimuthal moments, v n (p T ), shown in Fig. 4 which therefore differ dramatically between the two hadronization schemes. In particular, as a consequence of the almost complete flatness of the ∆φ distributions, both models present even harmonics for low-p T which are much smaller than the gluonic ones. To be more quantitative, in Tab. 1 we compare the magnitudes of the even azimuthal harmonics at the partonic and hadronic level. As one can see, the damping effect is almost total for p T < 1 GeV and persists up to p T 2 GeV.
Moreover, the prominent away-side peak present in the distributions of Fig. 3 causes large negative odd harmonics due to assumed local transverse momentum conservation.
It should also be noted that these v n arise solely from initial GB gluon bremsstrahlung and final hadronization effects. No final state interactions have been included in our simulations.

Hadronization mechanisms that damp gluon harmonics
At this point one needs to ask: what is the scale which determines the quenching of anisotropy for small values of the transverse momenta? Our picture has two natural scales: the typical energy of the hadronization process, λ 1 GeV, and the typical momentum exchanged, Q sat = 2 GeV. Which of the two plays the role of discriminant between the uniform and the anisotropic regions?
To answer this question we repeated the previous simulation using the independent fragmentation scheme but with a different value for the average momentum exchanged, Q sat = 4 GeV. The results are reported in Fig. 5. If compared to Fig. 2 one can appreciate that the final hadronic anisotropies now follow the gluonic ones much more closely than in the Q sat = 2 GeV case, as one might expect, for almost all values of the transverse momentum except for p T ≤ 1 GeV, where the two-pion distributions are again completely flattened out. This shows that the energy scale relevant in this decoherence effect is indeed the scale of hadronization, λ 1 GeV.
As a further check, one can study the average ratio between the pion and the gluon momenta both for the transverse and longitudinal components. Our simulation indicates that for the two models This again is a hint for the role of hadronization in damping the anisotropies at low p T . The reasons for this "quenching power" of the hadronization process can be found in essentially two features: (i) non-collinearity of the gluon radiation and (ii) isotropic decays of resonances. In particular, the fact that gluon fragmentation is not perfectly collinear is already manifest in the first order pQCD Altarelli-Parisi distribution, where k is again the gluon transverse momentum, z is the Bjorken variable and P(z) the splitting function. Although this distribution is strongly peaked around k 0 it also has very long tails, clearly showing that a totally collinear picture is a too naïve approximation (for an interesting discussion on the role of fragmentation functions on the azimuthal harmonics see [29]). Morever, in the LUND string model we also have the transverse fluctuations of the flux tube, whose typical scale is √ κ 0.45 GeV (see Eq. (4)), which provide another source of dilution of the initial anisotropy. Lastly, the intermediate steps between the initial gluons and the final pions are populated by resonances. These particles decay without any preferential direction and hence strongly contribute to the flattening of the final ∆φ distributions. In particular, this is the reason why the observed decoherence happens for p T 1 GeV, this being the typical mass of the most common resonances. To better illustrate this point we performed an overly simplified simulation whose description and results are shown in Tab. 2.
We learn from above that the strong damping of parton level v n at p T λ is not a consequence of the particular initial-state model chosen, but it is rather an intrinsic property of the hadronization mechanism. This is one of the striking results of our analysis.  Table 2: Output of two events composed by a qq pair and a gluon for both independent fragmentation and LUND model. The quark and antiquark are moving in opposite directions along z-axis, each of them with energy E q = 10 GeV. The gluon, instead, flies away along the x-axis (η = 0, φ = 0) with energy E g = 3 GeV. We only report those particles with pseudo-rapidity |η| < 1, i.e. close to the gluon in phase space. One immediately notices that even though the initial gluon has zero azimuthal angle, the fragmentation produces particles with φ 0. Every event always contains at least one resonance (particles in parenthesis) which then decays isotropically. Moreover, the fact that even the resonances (which are produced in the very first step of hadronization) have non-zero anzimuthal angle shows that the gluon fragmentation is non-collinear. The final pions are therefore widely spread in ∆φ. They also have p T < 1 GeV as expected.

Hadronization damping of triangular and higher gluon harmonics
As one can see from Fig. 4, in contrast to the LUND model, the independent fragmentation leads to very small odd harmonicsalbeit fluctuations due to small statistics effects. The reason is that for two color antennas the system is back-to-back symmetric, with just two qq pairs recoiling from each other. Since at high-p T no other correlations are introduced, the final pions inherit the symmetries of the initial partons. To check the consequences of hadronization on the purely geometrical odd moments -and to reproduce a slightly more realistic configuration -we performed a simulation involving three color antenna simulating three recoiling beam jets  (   Figure 6: Two-particle ∆φ distributions obtained from three initial qq pairs in a triangular configuration for both the initial-state gluons (blue) and the final pions (red). The hadronization has been performed using the independent fragmentation scheme.
conserving transverse momentum. Both the simulation and the analysis closely follow what we explained in Sec. 2 but now we implemented a third quark-antiquark pair such that the whole system conserves the total momentum, i.e. if Q 1 , Q 2 and Q 3 are the momenta exchanged by each of the three pairs then Q 1 + Q 2 + Q 3 = 0. In this case we only used the independent fragmentation scheme, which lacked the odd v n in the first place. The results for the ∆φ distributions and for the final Fourier moments are shown in Fig. 6 and 7.

Results
As one immediately notices, the initial gluons ∆φ distributions now have a contribution from odd Fourier components, v 2n+1 . Once again, in the p T λ region none of the initial information is preserved and the pions correlation functions are totally flattened. Even though we suffer from a lack of statistics, it is also clear that in the higher transverse momentum regime the final particles distributions keep following the gluonic ones. In terms of Fourier harmonics, even though the initial gluons have definitely sizeable v n , the final pions moments in the low-p T region are essentially zero (less than 10%) -see again Fig. 7. Hadronization is once again very effective in quenching the azimuthal asymmetry of the system. For p T λ the gluonic and hadronic v n have, instead, similar values.

Conclusions
Our analysis shows that the process of hadronization can lead to major distortions of the azimuthal harmonics up to p T = 3 GeV, which are of interest in the search for signatures of perfect fluidity in p+A systems. In particular, one important aspect seems to be general and modelindependent: at small values of the transverse momentum of the hadrons the complexity of the  hadronization itself -namely non-collinearity of fragmentation and isotropic resonance decays -greatly reduces the information contained in the initial-state partons, causing an almost total smearing of the final hadronic spectrum. As shown in Sec. 3, this is intrinsic to the hadronization itself and hence should be taken into account by every model with initial-state anisotropy at least with a theoretical error band estimated by testing several hadronization schemes.
Many models of initial-state anisotropy seem to claim that the only scale relevant for the generation of non-zero v n is the typical momentum exchanged, Q sat . Our analysis shows that this is only true above a second hadronization scale, λ 1 GeV. Hadronization scheme choice does matter, and those schemes that assume local collinear parton-hadron duality or pure collinear fragmentation without resonance production may over-estimate final-state hadron harmonics. It is interesting to note that assuming that fragmentation is perfectly collinear and described by any fragmentation function f (z) should actually lead to an enhancement of the initial gluon anisotropy since low-p T pions come from gluons with higher transverse momentum that naturally have greater azimuthal asymmetry. This is in striking contrast with the results of our JETSET simulations, where collinearity of hadronization is broken.
A second conclusion is that the two-particle correlations in the p T λ 1 GeV region can strongly depend on the chosen hadronization model. Specifically, while in the independent fragmentation scheme the initial parton anisotropy almost completely transmits anisotropy to the final hadrons, in the LUND string model a new source of correlation due to transverse momentum conservation is introduced, leading to a large away-side peak in the pion ∆φ distributions that is not due to back-to-back mini jets, which are not taken into account in the present simulations. If these features also survive to a system more complicated than ours, then this leads to two important remarks: first of all, whenever performing a Monte Carlo simulation to study multi-particle correlations one should always pay particular attention to the degree of model-dependence of the simulation itself since this might introduce important biases on the final conclusions. Secondly, if one assumes the description of hadronization in terms of a QCD color string as a fairly accurate one, then we might have an unexpected source of non-flow away side correlation for particles with p T > 1 GeV which would be difficult to deconvolute from pQCD di-jets by any experimental analysis.
We emphasize again that the properties of the final two-hadron correlations computed in this study, for the most part, are due neither to a particular initial-state mechanism nor to any collective transport or hydrodynamic effects: they are genuine consequences of the hadronization process that are scheme dependent. There is unfortunately no guarantee of universality of hadronization scheme.
Lastly, we also notice how, considering the observed mass dependence of final-state correlations [30], it would be interesting to repeate the previous analysis for different flavors and study how hadronization effects are affected by the particle masses. We leave this analysis to a future study.