Particle Production in pA Collisions at Mid-Rapidity in the Color Glass Condensate

: Particle correlations in small collisions systems, like proton–nucleus, lie at the core of the discussion about whether quark–gluon plasma is produced in small systems. Both initial and final state explanations have been essayed to describe such correlations. In this work, we focus on the initial state explanations provided by the quantum effects in the initial wave function of the incoming hadrons, in the framework of the Color Glass Condensate effective theory. We describe the formalism indicating the different inputs required for phenomenological applications. We compare the results from two different models, finding that the results for azimuthal harmonics agree qualitatively, but show quantitative differences, particularly at transverse momenta above the saturation scale.

An accepted explanation of the origin of the ridge observed at the RHIC in heavy ion collisions (HICs) is the collective flow due to strong final state interactions, usually described in the framework of relativistic viscous hydrodynamics; see the reviews [15,16].The azimuthal correlations observed in small collision size systems (both pp and pA) are also successfully described by relativistic viscous hydrodynamics [16,17] (for weak coupling final state explanations, see, e.g., [18] and the references therein).However, the hydrodynamic explanation of azimuthal asymmetries in small systems looks tenuous.This has triggered a many efforts to understand whether the final state particles carry the imprints of the partonic correlations that exist in the initial state.Several mechanisms have been suggested to explain the ridge phenomenon within the Color Glass Condensate (CGC) [19,20] framework (see [21] for a review).
One of the most-successful initial state explanations of the azimuthal correlations is based on the local anisotropy of the target fields [22][23][24], whose underlying physics mechanism can be explained as follows.In this model, the hadronic target is assumed to be composed of domains of oriented chromoelectric fields in the transverse plane.The size of these domains is of the order of 1/Q s , where Q s is the saturation momentum of the target.If the particles in the incoming nucleus are correlated, the transverse distance between these two particles is much smaller than 1/Q s .Therefore, the incoming particles scatter off the same domain in the target, and they acquire a common final momentum, which reflects the initial state correlations to the final state.Numerical studies based on this model were performed in [25,26].
Another mechanism that successfully describes the azimuthal correlations from the initial state point of view is known as the "glasma graph approach" [27][28][29].Even though this approach to two-particle correlations is shown to be very successful in describing many features of the data [30][31][32][33][34][35], the physics behind this approach is not clear.This issue was addressed in [36], where it was shown that a genuine quantum effect, namely the Bose enhancement of the gluons in the projectile wave function, leads to final state correlations within the glasma graph approach.In [37][38][39], it was also shown that another physical effect, known as the Hanbury-Brown-Twiss (HBT) correlations between the produced gluons, is present in the glasma graph approach.The correlations between quarks have also been studied, and it was shown that quarks in the initial state experience Pauli blocking due to their fermionic nature [40], which leads to (short range in rapidity) anti-correlations between hadrons originating from quarks.
The standard glasma graph approach is valid for dilute-dilute collisions.The extension of this approach to the dilute-dense situation applicable to pA collisions was studied in [41][42][43][44] by taking into account multiple scattering effects for double-and triple-inclusive gluon production 1 .A similar study was also performed for the double-inclusive quark [47] and heavy quarkonia [48,49] production in pA collisions.The extension of the glasma graph approach that includes the multiple scattering effects was also used in [50] to study the effects of the spatial eccentricity of the projectile shape on the second flow harmonics, in [51] to study the effects of fluctuations in the multiplicity of the produced particles, and in [52] to study the correlations between the total multiplicity and second flow harmonic, as well as the correlations between the mean transverse momentum and second flow harmonic.Finally, a specific model was used within the extended glasma graph approach to study three and four gluon production in [53].
In this paper, we focus on the computation of double-inclusive gluon production beyond the glasma graph approximation in pA collisions.In Section 2, we discuss the different contributions to correlated production, identifying the Bose enhancement and the HBT correlations following the work performed in [42].In Section 3, we discuss a specific model to compute double-inclusive gluon production following the work performed in [53].In Section 4, we study numerically the two-particle correlations and compare the results of the model described in Section 3 with a more-rigorous, but numerically involved implementation of the Lipatov vertices.Finally, in Section 5, we present a concise summary and discussion of the double-inclusive gluon production and correlations in pA collisions.

Double-Inclusive Gluon Production in Dilute-Dense Scattering at Mid-Rapidity
The two-gluon spectrum with transverse momenta k 1 , k 2 in pA collisions at midrapidity was computed in [54,55], and it reads 2 Here, ρ a (x) is the color charge density of the dilute projectile with a = 1, . . ., N 2 c − 1 the color index and U(x) is the adjoint Wilson line that accounts for the scattering of a gluon at transverse position x, which is defined as where A − (z + , x) is the background field of the target.⟨• • • ⟩ P(T) corresponds to projectile (target) averaging, and the Weizsacker-Williams field is given by The projectile averaging of the color charge densities is performed by adopting the generalized McLerran-Venugopalan (MV) [56,57] model with the Gaussian weight functional, which amounts to factorizing the color charge densities into a product of all possible Wick contractions.Moreover, the average of two color charge densities is given by the following general form 3 : Then, the projectile averaging of the four color charge densities that appears in Equation ( 1) reads On the other hand, the dipole and quadrupole operators are defined in the standard way as with the corresponding Fourier transforms defined as Q(x 1 , x 2 , x 3 , x 4 ) = q 1 q 2 q 3 ,q 4 e −iq 1 •x 1 +iq 2 •x 2 −iq 3 •x 3 +iq 4 •x 4 Q(q 1 , q 2 , q 3 , q 4 ) .( 9) Then, using Equation (5) for the projectile averaging and the definitions of the dipole and quadrupole operators given above, the double-inclusive gluon spectrum (Equation (1)) can be written as with where L i (k, q) is the Lipatov vertex defined as The contributions given in Equations ( 11)-( 13) correspond to the three different contractions given on the right-hand side of Equation ( 5) for the projectile averaging.The next step in the analysis of the double-inclusive gluon production is to understand the target averaging.For this purpose, we adopted the so-called area-enhancement model (which was originally introduced in Refs.[58,59] and used in Refs.[41,42]), which can be explained as follows.The computation of the double-inclusive gluon spectrum requires integration over four transverse coordinates of the Wilson line structure that appears as either a quadrupole or double-dipole operator in the spectrum.The leading contribution to the spectrum is expected to arise when as many points as possible are away from each other in the transverse plane.This configuration would maximize the result of the transverse integration, and any other configuration where the points are close to each other would be suppressed by the transverse area of the projectile.On the other hand, in order for the scattering matrix to be non-vanishing, the scattering objects must be color singlet operators with a small transverse size of the order of 1/Q s .Thus, the configuration that maximizes the spectrum would correspond to the case where four points are combined into two color singlet pairs and the distance between the pairs is large.Within this model, the target-averaged quadrupole and double-dipole operators can be approximated as follows: where, for the simplicity of the notation, we introduced d(x, y) ≡ ⟨D(x, y)⟩ T .Finally, assuming translational invariance of the target-averaged dipole operators: the double-inclusive gluon spectrum given in Equation ( 10), within the area-enhancement model, can be organized as with The final result for the double-inclusive gluon spectrum in pA collisions at midrapidity within the area-enhancement model is given in Equation (18).Now, we can discuss the physical origin of the various terms.First of all, it is straightforward to realize that the spectrum is symmetric under (k 1 , k 2 ) → (k 1 , −k 2 ), which is known as the "accidental symmetry of the CGC".The effects of this symmetry have been widely discussed in the context of the particle correlations in the literature, and it is argued that it can be broken by (i) adopting a nonlinear Gaussian approximation for the double-dipole operator [60], including density corrections to the dilute projectile [61][62][63] or (ii) including the subeikonal corrections [64][65][66][67] in the computation of the double-inclusive gluon spectrum.
In order to understand the underlying physics mechanism behind each term in the double-inclusive gluon spectrum, let us first elaborate on the function µ 2 (k, p), which is defined as the correlation function of two color charge densities of the projectile and can have the following general form: Here, the function T can be understood as the transverse-momentum-dependent distribution of the valence charges, while the function F can be interpreted as a soft form factor, which is maximal when its argument is zero and vanishes when its argument is away from zero.The factor R that appears in the argument of the soft form factor in Equation ( 22) is the radius of the projectile.Finally, it is important to remind that (k 1 − q 1 ) and (k 2 − q 2 ) are the momenta of the two gluons in the projectile, k 1 and k 2 are the momenta of the produced gluons in the final state, and q 1 and q 2 are the respective momentum transfers to the produced gluons from the target during the interaction.Using the properties of the soft form factor F, each contribution in the double-inclusive gluon spectrum can be interpreted as follows: • I 0 term: This contribution corresponds to the uncorrelated production, which is simply the square of the single inclusive spectrum.• I 1 -first term: The first term in Equation ( 20) is proportional to (23) where the soft form factor is peaked when the momenta of the two gluons in the projectile wave function are very close to each other.Therefore, this term contributes to the Bose enhancement of the gluons in the projectile.
• I 1 -second term: The second term in Equation ( 20) is proportional to where the soft form factor is peaked when the momenta of the two produced gluons are close to each other.Therefore, this term is the HBT contribution.• I 2 -first term: The first term in Equation ( 21) is proportional to where the soft form factor is peaked when the momenta of the gluons in the target that are transferred to the projectile are close to each other.Thus, this term contributes to the Bose enhancement of the gluons in the target.• I 2 -second term: The second term in Equation ( 21) is proportional to where the soft factor is peaked when the momenta of the gluons in the projectile are close to each other.Therefore, this is a contribution of the Bose enhancement of the gluons in the projectile.However, compared to the similar contribution in I 1 , these terms are suppressed by 1/(N 2 c − 1).Thus, the contribution to the Bose enhancement of the projectile gluons from the second term in I 2 can be interpreted as N c -corrections to the ones observed in I 1 .

A Specific Model for Double-Inclusive Gluon Production in Dilute-Dense Scattering
In this section, we review a specific approach to particle production in pA collisions that was originally introduced in Ref. [53] to use the computations of the three-and fourparticle production.In this approach, a diagrammatic method is introduced that uses the symmetries between different terms in the computation, which significantly reduces the number of terms that need to be calculated.Here, we will not discuss the details of this diagrammatic approach and refer the interested reader to Ref. [53].In the rest of this section, we focus on the two-gluon production in this approach and briefly discuss the derivation and the results.
The analog expression of Equation ( 27) for the double-inclusive spectrum in this approach is written as Here, averaging over the color charge densities in the projectile wave function is performed in the same manner as in Equation ( 4), but with a different normalization, and it is performed in momentum space, which simply reads In the standard MV model, the profile function of the projectile distribution μ2 is chosen to be a Dirac delta function.For a more-realistic analysis, we use a Gaussian distribution: which is peaked around k + q = 0, with B p the gluonic transverse area of the projectile.The projectile averaging of the four charge densities appearing in Equation ( 27) is factorized into product of the possible Wick contractions as in Equation ( 5).On the other hand, in Equation ( 27), target averaging is encoded in the effective scattering amplitudes, which are defined as with L i (k, q) being the Lipatov vertex defined in Equation ( 14) and U ab (y) the standard Wilson line in the adjoint representation given in Equation ( 2).On the target side, the multiple interactions result in dipole (Equation ( 6)) and quadrupole (Equation ( 7)) operators, as is usual for the double-inclusive gluon production.Target averagings of these operators are performed again within the area-enhancement model, i.e., quadrupole and doubledipole operators are factorized into the target-averaged single-dipole operators, as given in Equations ( 15) and ( 16), respectively.Moreover, in this computation, we adopted the well-known Golec-Biernat-Wüstoff (GBW) model [68,69] for the target-averaged dipole operators, which in momentum space reads The Lipatov vertices appear as a product of two, for each produced gluon in the spectrum: which suffer from infrared divergences.Usually, these divergences are regulated by introducing an infrared cutoff in the integrations over momenta.Here, we followed a different approach and assumed a "Gaussian regularization" in the collinear limit, where (k − q) 2 ∼ 0. This amounts to using the following expression for the product of two Lipatov vertices in order to regulate the infrared divergences: where ξ 2 is a parameter with dimensions of the momentum squared.
A comment is in order here.The main problem in using a Gaussian regularization for the product of two Lipatov vertices given in Equation ( 33) is that it only depends on the momentum of the projectile gluons, k i − q i , and not on the momentum of the produced gluons, k i .Therefore, the contribution included in Equation ( 33) is restricted to that when the gluon is emitted from the projectile color charge density and, then, interacts with the target.Thus, it is a missing part of the physics included in the standard Lipatov vertices.In this approach, the final momentum is acquired by the interaction with the target, which is suitable for the collinear limit of the projectile.In principle, in this limit, the "hybrid factorization" [70,71] is employed, and it corresponds to the forward production of partons near the fragmentation region.Even though our approach is more suitable for central production, with the approximation introduced in Equation ( 33), its validity is reduced to the forward limit, but not yet near the proton fragmentation region.In this region, the projectile partons are defined in terms of the Wigner functions [34,35,47,60,72].However, the standard Wigner function approach does not account for quantum correlation effects such as the Bose enhancement and HBT correlations discussed in the previous section, while our approach takes into account these quantum effects (see [53] for a detailed derivation).Thus, for two partons in the projectile, the joint Wigner function that we use reads Finally, by using the MV model for computing the projectile averaging via Equations ( 28) and ( 29), adopting the area enhancement argument given in Equations ( 15) and ( 16), using the GBW model (Equation ( 31)) for the dipole operator for computing the target averages, and using a Gaussian form for the product of two Lipatov vertices (Equation ( 33)) that mimics the Wigner function approach, the double-inclusive gluon spectrum can be organized as follows: with the classical contribution with the leading N c quantum correlations: and the subleading N c quantum correlations:

Two-Particle Azimuthal Harmonics
In this section, we proceed to undertake the computation of the azimuthal harmonics derived from the two-particle distribution given in Equation ( 1) 4 .To effectively investigate the azimuthal harmonics inherent in two-particle correlations, we adopted the cumulant method 5 .The primary objective of this method is to diminish the influence of what is commonly referred to as the "non-flow" correlation.This term encompasses contributions to the correlation function originating from processes divergent from genuine collective flow, such as resonance decays or jet correlations.The use of the cumulant method is instrumental in refining the definition of azimuthal harmonics, thereby providing a moreaccurate characterization of the true collective flow dynamics.Within this framework, the 2-particle azimuthal harmonics of order n read as follows: where κ n {2} is the nth harmonic distribution of the 2-particle production: In the model introduced in Section 3, specifically through the Gaussian regularization of the Lipatov vertices in Equation (33), we can perform an analytical computation of two-particle azimuthal harmonics, as presented in Ref. [53].In Figure 1, we plot the relationship between v n {2} (n = 2, 4) and the dimensionless variables B p Q 2 s and ξ 2 /Q 2 s , keeping N c = 3 constant.Notably, at a fixed Q s , the behavior of even azimuthal harmonics undergoes discernible changes as we vary the proton area and ξ.A noteworthy observation is the rapid growth of even azimuthal harmonics when both B p and ξ approach zero.Conversely, a gradual decrease is observed as these parameters attain larger values.
This reduction in even azimuthal harmonics with increasing B p aligns with our expectations in the color domain model of particle correlation.The rationale behind this trend lies in the influence of the number of domains, denoted by B p Q 2 s .As this quantity increases, the likelihood of two gluons undergoing scattering within the same domain diminishes, resulting in an overall decrease in correlation.Likewise, the azimuthal harmonics v n can be expressed as a function of the transverse momentum.To achieve this, integration is performed over one of the particle transverse momenta in Equation ( 42), holding the other constant.Consequently, the computation and definition of the differential azimuthal harmonics are given as follows: Analogous to the integrated azimuthal harmonics, the use of the Gaussian regularization outlined in Equation ( 33) facilitates an analytical computation of the differential harmonics.The outcomes of v 2n {2}(p ⊥ ) for n = 1, 2, 3 are depicted in Figure 2, where B p = 6 GeV −2 , ξ = Q s /2, Q 2 s = 2 GeV 2 , and N c = 3.While our primary focus does not involve a direct comparison with experimental data, it is noteworthy that the obtained values qualitatively agree with empirical observations.In fact, the experimental values for v 2 {2} lie in the range 10-20% for p ⊥ ∼ 1-2 GeV/c [6][7][8][9]12,13].It is essential to highlight that, given the Gaussian forms employed, our results may not be deemed reliable for p ⊥ significantly exceeding Q s .Finally, we undertake a comparative analysis of the azimuthal harmonics, computed using the Gaussian regularization applied to the Lipatov vertices as detailed in Equation (33), with those derived from the same computation using the actual expression presented in Equation (32), albeit subject to regularization: where the parameter µ 2 g serves as an effective mass introduced to address and regularize the infrared divergences 6 .
To conduct the comparison, we fixed the values for the parameters: µ g = 0.5 GeV, Q 2 s = 2 GeV 2 , ξ 2 = 0.5 GeV 2 , B p = 6 GeV −2 , and N c = 3.We start by a comparative analysis of the integrated momentum distribution, defined as follows: In Figure 3, we compare the distributions obtained using the regularization schemes specified in Equation ( 33), denoted as Gaussian, and Equation ( 44), denoted as Lorentzian.Notably, both results exhibit comparable magnitudes; however, the Lorentzian regularization manifests a more-restrained tail, aligning with our expectations (see the comments on the lack of dependence on the transverse momentum of the produced gluons below Equation ( 33)).This observation suggests that, with Gaussian regularization, there is an overestimation of high p ⊥ particles, indicating the unreliability in this regime.(45), is illustrated as a function of the transverse momentum.The blue line corresponds to the spectrum employing Gaussian regularization, while the red line corresponds to the spectrum using Lorentzian regularization.
In Figure 4, the two-particle integrated azimuthal harmonic is depicted as a function of B p employing both Gaussian and Lorentzian regularizations.Notably, as the projectile area increases, both results exhibit a decrease, yet the Gaussian regularization yields harmonics approximately four-times larger.A similar trend is observed in Figure 5, where the differential harmonics are plotted against the transverse momentum.While both results follow a comparable trend, the Gaussian regularization generates significantly higher harmonics.
While both approaches, Gaussian and Lorentzian, give results that lie in the ballpark of the experimental data as mentioned before, our study shows that large uncertainties remain.Besides, we should take into account other effects like fragmentation, the contribution from the quark channel, and the possibility of final state effects.Also, non-eikonal contributions may be sizable at RHIC energies.All these considerations prevent us from a quantitative comparison to experimental data.

Conclusions
In this manuscript, we described the approach to two-particle correlations in the Color Glass Condensate effective theory.After describing the formalism and the physical mechanisms at work in two-particle inclusive distributions, we went through the different model-dependent ingredients required for phenomenological applications.We presented a Gaussian model developed in [53], based on the use of Wigner functions, that allows analytic computations.This model oversimplifies the treatment of the large transverse momentum tails in particle production and should be valid for transverse momenta of the order of the saturation scale.
Then, we numerically computed the azimuthal harmonics v n=2,4,6 and studied their dependencies on different parameters in the model.Finally, we compared the results of this model with a more-numerically involved implementation, which keeps the perturbative form of the Lipatov vertex, albeit with a regularization in the IR.We showed that this implementation results in smaller azimuthal harmonics, particularly at the transverse momentum larger than the saturation scale.We concluded that the uncertainties resulting from the details of the modeling required for phenomenological applications may be significant and have to be considered for a meaningful comparison with experimental data.

v 4 2 Figure 1 .
Figure 1.The relationship between the even 2-particle azimuthal harmonics, represented as v n {2}, is examined in terms of its dependence on ξ 2 /Q 2 s (left panel) and B p Q 2 s (right panel).The left panel presents the variation of azimuthal harmonics while keeping B p Q 2 s = 12, while the right panel explores the impact of changes in B p Q 2 s with a fixed value of ξ 2 /Q 2 s = 1/4.

Figure 3 .
Figure 3.The single-particle momentum distribution, as expressed in Equation(45), is illustrated as a

v 2 2 Figure 4 .Figure 5 .
Figure 4.The 2-particle azimuthal harmonic is presented as a function of B p .The blue line corresponds to the spectrum employing Gaussian regularization, while the red line corresponds to the spectrum using Lorentzian regularization.

Author
Contributions: Conceptualization, P.A., T.A. and N.A.; methodology, P.A., T.A. and N.A.; investigation, P.A., T.A. and N.A.; writing-original draft preparation, P.A., T.A. and N.A.; writing-review and editing, P.A., T.A. and N.A.All authors have read and agreed to the published version of the manuscript.Funding: P.A. and N.A. have received financial support from Xunta de Galicia (Centro singular de investigación de Galicia accreditation 2019-2022, ref.ED421G-2019/05), by the European Union ERDF, by the "María de Maeztu" Units of Excellence program MDM2016-0692, and by the Spanish Research State Agency under project PID2020-119632GBI00.This work was performed in the framework of the European Research Council project ERC-2018-ADG-835105 YoctoLHC and the MSCA RISE 823947 "Heavy ion collisions: collectivity and precision in saturation physics" (HIEIC) and has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No. 824093.P.A. has been supported by Consellería de Cultura, Educación e Universidade of Xunta de Galicia under the grant ED481B-2022-050.