Forward di-hadron back-to-back correlations in $\boldsymbol{pA}$ collisions from rcBK evolution

We study the disappearance of the away-side peak of the di-hadron correlation function in p+A vs p+p collisions at forward rapidities, when the scaterring process presents a manifest dilute-dense asymmetry. We improve the state-of-the-art description of this phenomenon in the framework of the Color Glass Condensate (CGC), for hadrons produced nearly back-to-back. In that case, the gluon content of the saturated nuclear target can be described with transverse-momentum-dependent gluon distributions, whose small-$x$ evolution we calculate numerically by solving the Balitsky-Kovchegov equation with running coupling corrections. We first show that our formalism provides a good description of the disappearance of the away-side azimuthal correlations in d+Au collisions observed at BNL Relativistic Heavy Ion Collider (RHIC) energies. Then, we predict the away-side peak of upcoming p+Au data at $~\sqrt[]{s}=200$ GeV to be suppressed by about a factor 2 with respect to p+p collisions, and we propose to study the rapidity dependence of that suppression as a complementary strong evidence of gluon saturation in experimental data.


I. INTRODUCTION
Azimuthal correlations of particles in the final states of hadronic collisions serve as a powerful tool for experimental tests of the Color Glass Condensate (CGC) [1][2][3], the effective theory of protons and nuclei in the nonlinear regime of quantum chromodynamics.A special role in the phenomenology of the CGC is played by correlations of particles in p+A collisions probed in the region of fragmentation of the protons [4][5][6][7][8][9][10][11], where the rapidities of the correlated particles are large and positive (forward rapidity region).Such configurations are ideal for testing the CGC theory, because they induce a dilute-dense asymmetry in the problem: The projectile proton is probed at large values of Bjorken x, and is thus a dilute object, amenable to a description in terms of well-known parton distribution functions (PDFs); The nuclear target is instead seen as dense state of low x gluons, a regime in which the saturation of the gluon densities is manifest, so that the CGC description applies.This dilute-dense asymmetry hence minimizes our uncertainty in the knowledge of the projectile, and provides the cleanest possible environment for the study of phenomenological signatures of gluon saturation in the target.
In this paper, we deal with a salient prediction of the CGC theory: The disappearance of the away-side peak (∆φ = π) of the two-particle correlation function of dilute-dense collisions (i.e., forward p+A collisions).Following [12], let us provide an intuitive picture of this phenomenon.A valence parton interacting with a CGC (i.e., a large classical Yang-Mills background field) undergoes multiple scattering with low-x gluons, either before or after splitting into a pair of back-to-back partons, which eventually produce the jets or hadrons observed in the final state.The pair of partons is put on shell via the interaction with the target, and this occurs through a transverse momentum exchange of order of the saturation scale of the target, Q s , which is typically much larger than the transverse momentum of the parent valence parton.The back-to-back correlation of the final-state particles, which would not be affected by an interaction with a gluon of zero transverse momentum, is therefore altered, and this induces a depletion of the correlation function around ∆φ = π.Consequently, the away-side peak observed in p+A collisions is expected to be suppressed with respect to that of p+p collisions, because the target nuclei are denser and more saturated.
Experimentally, the validity of this picture is strongly supported by Relativistic Heavy-Ion Collider (RHIC) data, as both the STAR and the PHENIX collaborations reported a visible suppression of the away-side peak when comparing p+p collisions to central d+Au collisions [13,14].These data, though, suffer from large uncertainties.More accurate tests of the CGC prediction may nevertheless become possible with the advent of data from the recent 200 GeV p+Au run performed at RHIC.As we shall see, one of the goal of this paper is to provide predictions for the away-side peak in these collisions.
On the theory side, first calculations of forward two-particle production in p+A collisions within the CGC framework date back to more than 10 yrs ago [4,15].The cross section for the production of two particles is intrinsically difficult to evaluate, because it involves multi-point correlators of Wilson lines.Over the years, different levels of approximation have been employed to perform calculations and obtain predictions, as reviewed in Ref. [12].The simplest option is to disregard non-linear effects and recover the so-called k t −factorization (or high-energy factorization) framework [5,16]; the cross-section is then obtained from a single two-point correlator, but that approximation is not applicable in the away-side peak region.In Ref. [6], the multi-point correlators are evaluated using the so-called Gaussian approximation of the non-linear QCD evolution, however only the elastic contributions are kept, and it turns out that the neglected contributions are also sizable in the away-side peak region.In Ref. [8], the complete Gaussian expressions are used, however due to the complexity of the problem, only quark-initiated channels could be included.
A crucial step was the realization that the cross section simplifies dramatically if one considers the production of partons which are nearly back-to-back [17,18].In this limit, the dense part of the scattering (the nucleus) is characterized by transverse-momentum-dependent (TMD) gluon distributions whose small x evolution is easily affordable to numerical implementations, because the multi-point correlators of Wilson line [19] involve only two distinct transverse positions.This framework has been employed in a number of applications [7,10,11,[19][20][21][22], and was recently reviewed in [23] 1 .In the case of forward di-hadron production, it has only been used together with Golec-Biernat Wusthoff (GBW) type parameterizations for the gluon TMDs [7], which suffer from unphysical exponential tails at large gluon transverse momentum.The goal of our work is to improve on this by obtaining the TMD gluon distributions from numerical solutions of the QCD non-linear evolution.
In [19,20], the gluon TMDs and their small-x evolution were obtained from the full QCD evolution at leading logarithmic accuracy, i.e. from the Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) equation.Since the implementation of running-coupling corrections in this context has not been performed yet, we prefer to work within the Gaussian approximation of JIMWLK evolution and obtain the gluon TMDs from the Balitsky-Kovchegov (BK) equation with running coupling corrections (rcBK), because we expect running-coupling corrections to be much more important than corrections to the Gaussian approximation.In addition, the rcBK solutions are well constrained from deep inelastic scattering data [25], so that the final expression of the cross section of two-parton production turns out to be essentially free from tunable parameters.We shall derive this cross section, and convolute it with fragmentation functions to present state-of-the-art results on azimuthal correlations of di-hadrons in forward p+A and p+p collisions at RHIC energies.We both test our theory against existing data and make predictions for future back-to-back correlations of hadrons.
The paper is organized as follows.In Sec.II, we briefly review the theoretical formalism of nearly-back-to-back forward di-hadron production in p+A collisions in the CGC framework, and we present the fully differential cross section for the production of di-hadrons, specifically, two neutral pions.In the cross section we shall introduce the TMD gluon distributions which characterize the dense component of the scattering process.In Sec.III, we explain in detail how such quantities are obtained from rcBK evolution, and we show their behavior as function of the kinematic variables.Calculations of the away-side peak are eventually given in Sec.IV.The per-trigger-yield cross section is calculated as function of the relative azimuthal angle of the two hadrons.We first calculate it in d+Au and p+p collisions, and we compare our results to existing RHIC data.Then, we compute several predictions for the away-side peak in upcoming p+Au collisions at √ s = 200 GeV.In Sec.V, we make predictions for the evolution with rapidity of the suppression of the away-side peak, and we compare our results with calculations performed using an alternative implementation of the rcBK evolution for nuclei.Section VI is left for conclusive remarks.

II. COLOR GLASS CONDENSATE IN THE BACK-TO-BACK REGION: TMD FACTORIZATION
We study the production of pairs of hadrons in forward p+A collisions.We display such process in Fig. 1.Working in light-cone coordinates, (+, ⊥, −), the Feynman-x variables associated with the projectile parton moving along the + direction and with the target gluon coming from the − direction are, respectively, given by where the k i 's refer to the outgoing partons while the p i s refer to the final-state hadrons.We have introduced their transverse momenta, p 1t = z 1 k 1t and p 2t = z 2 k 2t , and rapidities, y 1 and y 2 , while √ s denotes the invariant mass of the scattering process.Equation (1) shows that when both y 1 and y 2 are large and positive, we probe a large-x parton in the projectile, and a small-x gluon inside the target nucleus.In experiments at RHIC, one can reach y ∼ 4, leading to x 1 ∼ 0.5, and x 2 ∼ 10 −3 : This realizes the anticipated dilute-dense asymmetry of forward particle production, which is essential for the applicability of our formalism.Now, following [18,19], we dub and we introduce the following variables If we stick to a limit in which the produced particles are back-to-back, i.e., their relative azimuthal angle, ∆φ, is close to π, then the total transverse momentum of the di-hadron pair is much smaller than the transverse momentum of the single hadrons, i.e., |k t | |P t | [18].
Following the exhaustive derivations of [19], the advantage of this limit is that it allows to write the cross section of the scattering process as an expansion in powers of 1/P t .Keeping only the leading order terms in this expansion, the dense component of the scattering is given by a combination of transverse momentum dependent gluon distributions (TMDs in short), which are CGC correlators of traces of Wilson lines.Summing over all production channels (qg → qg, gg → q q, gg → gg), the cross section for the production of two partons can be written in the following compact notation [19] where we note the manifest factorization of the cross section into a dilute component, characterized by collinear parton distribution functions f a/p (x 1 , µ 2 ), evaluated at a factorization scale µ 2 , and a channel-dependent dense component, characterized by hard factors [18], H(z, P t ), and the TMD gluon distributions, F (i) (x 2 , k t ), specified below.
To turn Eq. ( 4) into a tool enabling us to compute predictions for di-hadron production, we convolute it with fragmentation functions.Considering u quarks, d quarks and gluons in the projectile proton, and considering only their fragmentation into pions, and neglecting all terms which are suppressed by 1/N 2 c , the full expression of the cross section reads where we have denoted by D π 0 /a (z i , µ 2 ) the fragmentation of a parton a into a neutral pion at the factorization scale µ 2 , and the notation used for the distributions is the same as in [19].In Sec.IV we shall make use of Eq. ( 5) to compute azimuthal correlations of neutral pions at RHIC.Let us first describe, in the following section, the rcBK formalism developed for the small x 2 evolution of the TMD distributions which appear in the cross section.

III. EVOLUTION OF THE TMD GLUON DISTRIBUTIONS TOWARDS SMALL x
In order to complete our formulation of the cross section, Eq. ( 5), we discuss now the x 2 evolution of the TMD gluon distributions, F (i) ag (x 2 , k t ).The starting point is the evolution of the impact parameter (b) independent fundamentaldipole scattering amplitude, which we denote in a standard notation N F (x, r).As it is customarily done in the literature, we assume that the b dependence of N F factorizes, and that it does not mix with the evolution.The evolution equation of the dipole amplitude, known as the Balitsky-Kovchegov equation [26,27], supplemented with running coupling corrections (rcBK equation), reads (r i = |r i |) with r 2 ≡ r − r 1 , and where x 0 is some initial value for the evolution (usually chosen to be x 0 = 0.01).K run is the evolution kernel including running coupling corrections.Different prescriptions have been proposed in the literature for K run .As shown in [28], Balitsky's prescription minimizes the role of higher conformal corrections: The rcBK evolution is independent of whether the target is a proton or a nucleus.That is accounted for in the initial condition.We use the so-called McLerran-Venugopalan (MV) model: with Λ = 0.241 GeV, and where Q s0 denotes the saturation scale at the initial value x 0 .We use x 0 = 0.01 and Q 2 s0 = 0.2 GeV 2 for a proton target, which are known to provide a good description of single-inclusive forward hadron RHIC data [29].For a target nucleus, things are more uncertain, as we are interested only in central collisions, i.e., collisions at small impact parameter.Motivated by previous studies [6], we keep x 0 = 0.01, and we choose Q2 s0 = 0.6 GeV 2 , i.e., a factor 3 larger than the Q 2 s0 with a target proton.Now, the simplest gluon TMD distribution, F qg (x 2 , k t ), is related to the Fourier transform of the fundamental dipole amplitude, N F (x 2 , r), and is given by: where and with S ⊥ denoting the transverse area of the target.
In full generality, none of the other gluon TMDs can be obtained in such a straightforward manner, directly from N F , or its Fourier transform F .To move forward, we resort to a mean-field type approximation: we shall utilize the so-called Gaussian approximation of the CGC [4,[30][31][32][33][34][35].The essence of this approximation is to consider all the color charge correlations in the target to stay Gaussian throughout the evolution.This approximation, along with the large N c -limit, ensures the factorization of CGC expectation values into single-trace expectation values, and allows to calculate gg and F (2) gg from F [18]: We note that the difference which enters the cross-section (5), is related to the adjoint-dipole scattering amplitude N A in the same way that F (1) qg was related to the fundamental dipole scattering amplitude.Indeed, if we introduce where then one has F (1) gg = F adj .This identify is true in full generality, beyond the Gaussian and large-N c approximations (for which 1 ) used here, as was first noticed in [20].
Finally, the two remaining gluon TMDs need to be computed from the Weizsäcker-Williams (WW) gluon distribution [18], which we denote F W W , and which should be obtained from the quadrupole operator Tr [A(x)A(y)] x2 where A(x) = U † (x)∂ x U (x) with U denoting a Wilson line.Again, this quantity is in general not related to the solution of the BK equation, F (x 2 , k t ), but using the Gaussian approximation one can write (in the large N c limit) 2 : F gg F gg (1) -F gg F gg Y = 0 FIG. 2. This figure presents the x2 evolution of three TMDs appearing in the cross section of Eq. ( 5), for a target proton.In panel (a), the initial conditions at x2 = 0.01 are presented.We show F (1) gg − F (2) gg (dashed line), and F gg (dotdashed line).The vertical dotted lines represent the saturation scale at the given value of x2.In the figure, Y = ln 0.01/x2 .The plotted quantities do not include the factor S ⊥ /αs in Eq. ( 9), common to all the gluon TMDs.
which allows to calculate the remaining two gluon TMDs needed in the cross section as follows [18]: We have now expressed all the needed gluon TMDs in terms of F (x 2 , k t ), the solution of the BK equation.We show in Fig. 2 some of those gluon distributions for a target proton, as function of k t , and for two values of x 2 .We do not show F (2) gg explicitly, but rather the difference F (1) gg , which effectively plays a role in Eq. ( 5).The TMD distributions present three specific features, fully characterizing the dense component of our scattering.Starting from the region where k t 1 GeV, we note that all the curves approach the same asymptotic behavior, i.e., an inverse power law, precisely equal to k −2 t at x 2 = x 0 [panel (a)], with a smaller absolute slope after x 2 evolution [panel (b)].As k t becomes of order 1 GeV, the TMD distributions start to separate, and quickly change their slope at a specific k t , which corresponds approximately to the location of the maximum of F (1) gg .This is the value of the saturation scale, Q s , which we indicate in both panels with a vertical dotted line 3 .Note that the small-x 2 evolution has the effect of shifting the saturation scale to larger values.Eventually, below the value of Q s the distributions become flat, and saturation is manifest.The difference F (1) gg goes to zero at k t = 0, consistently with Eq. ( 14).

IV. THE AWAY-SIDE PEAK FROM rcBK EVOLUTION: RESULTS AND PREDICTIONS
We can eventually employ the theoretical formalism introduced in the previous sections to compute azimuthal correlations of two hadrons in p+p and p+A collisions at √ s = 200 GeV.We integrate Eq. ( 5) over the momenta and rapidities of the produced hadrons, and study the behavior of the cross section as function of relative azimuthal angle, ∆φ.In the notation of [6], the observable we want to calculate is 3 More specifically, the maximum of F adj corresponds to the adjoint saturation scale, which is 1.5 times bigger than Qs, the fundamental saturation scale which corresponds to the maximum of F qg .This explains why the vertical line in Fig. 2(a) does not correspond to The experimentally measured quantity is not directly given by Eq. (20).Experimentalists normalize N pair (∆φ) with the total number of hadrons that trigger the correlations, i.e., in which we have introduced the cross section for single hadron production [29] dσ pA→π 0 +X dy dp where F and F are computed from the rcBK evolution equation as explained in the previous section.The final observable is dubbed coincidence probability by the STAR collaboration, and is given by Before showing our results, let us list all the details about the quantities needed in the calculation of CP (∆φ): • The parton distribution functions (PDFs) describing the projectile are taken from the NLO MSTW2008 fits [36]; • The fragmentation functions (FFs) used are the recent DSS14 NLO sets [37]; • The strong coupling constant appearing in Eq. ( 5) is calculated at NLO, and is given by the following expression where we take N f = 4, and Λ = 197 MeV.For µ 2 , we use the same scale employed in the PDFs and in the FFs (see item below); • The PDFs, the FFs, and α s are computed at the scale µ 2 = p 2 t1 , i.e., at the transverse momentum of the leading hadron.

A. Comparison with Run-8 d+Au RHIC data
Saturation effects are expected to yield a larger CP (∆φ) in pp collisions than in pA collisions, when ∆φ is in the vicinity of π.STAR data on CP (∆φ) [13] for neutral pion correlations are shown as symbols in Fig. 3. Data present a visible suppression of the correlation in d+Au collisions, suggesting that saturation effects may be effectively at play.The outcome of integrating Eq. ( 5) over the STAR kinematics in both p+p and d+Au collisions4 is shown as shaded bands, in the nearby of ∆φ = π.The shaded bands represent the uncertainty in the choice of the factorization scale, µ 2 , in Eq. ( 5).The upper limit of the bands is obtained with µ 2 = p 2 t1 , and the shaded area is obtained by taking a scale larger by 50%. 5 In the following, we shall therefore always employ µ 2 = p 2 t1 , which should provide the best agreement with data.
Figure 3 shows that the suppression of the away-side peak provided by our calculation is in agreement with the data, although robust conclusions are impossible to draw due to the large uncertainties in d+Au collisions.We also notice that our calculation reasonably captures the magnitude of CP (∆φ) at the away-side peak of p+p collisions.What we fail in reproducing, though, is the width of the measured correlation in p+p, which appears to be broader than our result.This has a simple explanation: in our calculation we are not supplementing the cross section with Sudakov factors, i.e., we do not take into account the radiation of soft gluons in both the initial and the final state, which would naturally provide a broadening of the correlation function.This was done in Ref. [38], using GBW-type parametrization for the gluon TMDs, but it remains to be done with the rcBK gluon TMDs.An attempt of such Sudakov resummation in the case of di-jet production was made within the Kutak-Sapeta (KS) approach [39] (see Sec. V), and the results are promising, in the sense that one observes a clear broadening around ∆φ = π when Sudakov resummation is included.We finally note that the correlation function shown in Fig. 3 is somewhat less flat than the one obtained in [29].By comparison, our formalism is valid in a narrower window near ∆φ = π, but it is more accurate there.

B. Predictions for p+Au collisions
In Fig. 4 we present predictions for the away-side peak of neutral pions in p+p and p+Au collisions at √ s = 200 GeV.This is achieved by integrating Eq. ( 5) over the kinematic cuts used by the STAR collaboration in their new analysis.We predict that the away-side peak is suppressed in p+Au by a factor close to 2. We find this conclusion to be rather independent of the p t window chosen for the measurement.

V. RAPIDITY DEPENDENCE OF THE SUPPRESSION AND COMPARISON WITH THE KUTAK-SAPETA APPROACH
A generic prediction of the CGC framework is that any effect due to gluon saturation should become less visible if we move towards more central rapidities, i.e., in our case, if we reduce the dilute-dense asymmetry by probing larger values of x in the nuclei.Consequently, the suppression of the away-side peak in p+A collisions relative to p+p should essentially fade away if we correlate particles in more central rapidity intervals.It is important to stress that the dependence on rapidity is a very specific feature of the saturation framework, which is not predicted by typical competing effects, e.g., conservation of total transverse momentum [40], or other energy-momentum conservation corrections which are relevant in the proximity of x 1 → 1 [41,42].Another competing description is that reported by Kang et al. [43], who manage to describe the suppression of the away-side peak without resorting to a CGC description, but solely from (cold) nuclear transverse-momentum broadening effects.Such models do not predict a specific dependence on the rapidity, so that the CGC interpretation would be strongly favored if such dependence is observed in data.Let us stress that the away-side peak in different rapidity intervals could be easily measured at the STAR or at the LHCb detectors, which present wide rapidity coverages.
Let us show, then, what our formalism predicts for the rapidity dependence of the suppression of the away-side peak.For reasons which will appear clear in the following discussion, it is very instructive to perform calculations and show results using both our rcBK formalism and the alternative Kutak-Sapeta (KS) approach [16], which we briefly review below.
In the KS approach, the momentum space version of the BK equation is used (written below for F p = πF qg , for a target proton): This way of writing the BK equation is convenient as it allows to include relatively easily some higher-order corrections, and in particular running-coupling corrections [44].To write down the non-linear term of Eq. (25) (last line in the equation) for the impact-parameter-integrated gluon distribution, it is assumed that integration over impact parameter yields d 2 b = πR 2 , where R is the radius of the target proton.The evolution of the gluon TMD in the case of a nucleus, F A , is then obtained through the following formal substitution in Eq. ( 25), In the above equation, R A is the nuclear radius, A is the mass number (A = 208 for Pb), and c is a parameter that is supposed to vary between 0.5 and 1, to assess the uncertainty related to the nonlinear term.The density F A obtained from Eq. ( 25) with the substitution above is the nuclear gluon density normalized to the number of nucleons in the nuclei.
The KS evolution in Eq. ( 25) is A-dependent through the non-linear term (it has to be so, since F A is an impact parameter integrated distribution), but the prescription for the initial condition is to choose the same in the nuclear case as in the proton case, i.e., F A (x 0 , k 2 ) = F p (x 0 , k 2 ).This is the major difference with respect to the approach presented in Sec.III, where an A-dependent initial condition and an A-independent evolution were used.
Figure 5 shows an illustration of the effect due to this difference between the rcBK and the KS approaches, which are both based on the same small-x evolution.In the figure we show F (1) gg for a target nucleus divided by the same quantity for a target proton 6 , for different values of rapidity Y , which is defined as Y = ln 0.01/x 2 .On the left, the rcBK distributions predict the same amount of suppression at each value of Y in the fully saturated region, k t ∼ 0, because the small-x evolution is A-independent.This is not the case in panel (b), where the ratio in the KS scheme is equal to unity for Y = 0 (not shown), and the difference in the evolution of the nucleus with respect to the proton is manifest already at k t = 0. Note that the plot is drawn for c = 0.5, but we stress that the qualitative picture is essentially independent of the choice of this constant.
This difference shown in Fig. 5 has a non-negligible impact on the rapidity dependence of the suppression of the away-side peak, which is the subject under study in this section.To show this, we calculate the ratio CP (∆φ) in p+Au over the same quantity in p+p using both the standard rcBK approach and the KS alternative proposal7 , and we look at its dependence with rapidity.We keep the old STAR kinematics of Fig. 3 for the p t of the produced hadrons, and we compute CP (∆φ) pA /CP (∆φ) pp , with an obvious meaning of the notation, around ∆φ = π in different intervals of rapidity.Results are shown in Fig. 6.We find that both schemes provide a hierarchy as function of rapidity expected in the saturation framework: To more central rapidities correspond larger values of the ratio around ∆φ = π, i.e., less suppression at the away-side peak.We stress that this is a peculiar feature of the saturation framework, and we strongly encourage measurements of this ratio in different rapidity intervals, which could provide, arguably, the strongest possible evidence in favor of the saturation picture.In addition, we expect such quantity to be almost unaffected by the uncertainties on the factorization scale (which turned out to be quite large in Fig. 3), as they are likely to cancel in the ratio.
Besides confirming the generic prediction of the CGC framework, precise measurements in p+Au collisions might as well shed light on the very validity of the approaches taken for the small-x evolution of the dense targets.In Fig. 6 we observe two notable differences between rcBK and KS.First, the dependence on rapidity at ∆φ = π is about twice stronger in the KS approach [panel (b)]: This results from having a small-x evolution at low k t (Fig. 5).Second, the rcBK case presents ratios which grow towards unity much faster as we move away from the back-to-back region.Specifically, the ratio at ∆φ = 3 is larger by 15% in the rcBK scheme.Such visible differences are expected to be sizable in the upcoming data, and would help improve significantly our understanding of the evolution equations of QCD in the nonlinear small-x regime.

VI. CONCLUSION
We have calculated the production of back-to-back pions in p+A and p+p collisions at RHIC energies, using the state-of-the-art CGC framework, i.e., the cross section reported in Eq. ( 5).We have developed a novel approach for the small-x evolution of the TMD gluon distributions F (i) ag , in which they are obtained from the BK evolution with an evolution kernel that includes running coupling corrections.The evolution is identical for proton and nuclear targets, the only difference being the value of Q 2 s at the initial condition.The validity of our framework is confirmed by the good agreement observed between the available data and our results in Fig. 3.
We thus derived genuine predictions of the CGC theory.The away-side peak in upcoming p+Au data is suppressed by about a factor 2 with respect to p+p collisions (Fig. 4), and this suppression tends to disappear as we reduce the dilute-dense asymmetry of the problem (Fig. 6).We stress, once more, that the combination of these two effects is a much stronger probe of gluon saturation than the suppression of the away-side peak alone.We have further compared the expectation of our framework to those of another state-of-the-art rcBK implementation, namely, the KS approach.Using the observable proposed in Fig. 6, p+Au data will potentially allow us to make a data-driven distinction between these two schemes of small-x evolution.
Before concluding, we stress that our calculation lacks an important ingredient: The inclusion of the soft gluon resummation, i.e., of Sudakov factors attached to the cross section which could potentially solve our problem of a too narrow correlation peak around ∆φ = π (Fig. 3).This improvement of our formalism will be presented in an upcoming publication.

FIG. 1 .
FIG.1.The pA → π 0 π 0 X process.See the text for details about the displayed quantities.

FIG. 3 .
FIG.3.The figure shows STAR data on azimuthal π 0 correlations at forward rapidity, in p+p collisions (circles) and central d+Au collisions (triangles) at √ s = 200 GeV.To remove fake two particle correlations which are essentially due to pileup effects, an arbitrary offset is added to push the STAR measurements close to 0 at the minimum of the correlation functions.Calculations of CP (∆φ) in our TMD+rcBK framework are shown as shaded bands.Light-shaded band: p+p collisions.Dark-shaded band: d+Au collisions.The meaning of the shaded bands is discussed in the text.

FIG. 4 .
FIG.4.In this figure we show predictions for azimuthal correlation of forward neutral pions in p+p (dashed line) and p+Au (dotted line) collisions at √ s = 200 GeV.Different panels correspond to different pt cuts applied to the cross section.

FIG. 5 .
FIG. 5.The figure shows F (1) gg for a target nucleus divided by the same quantity for a target proton, as function of kt.Results are shown within two different evolution schemes, namely rcBK [panel (a)] and KS approximation [panel(b)].The ratio is taken at different values of x2, indicated with different line styles.In the figure, Y = ln 0.01/x2 .

FIG. 6 .
FIG. 6.The figure shows the ratio CP (∆φ)pA/CP (∆φ)pp around ∆φ = π.Different line styles represent different rapidity intervals.Panel (a) shows results with gluon TMDs obtained as described in section III.In panel (b) the TMDs are obtained using the KS scheme, with c = 0.5.