Single spin asymmetry in forward $pA$ collisions: Phenomenology at RHIC

We confront the theoretical result of single spin asymmetry (SSA) $A_N$ in forward $pA$ collisions $p^\uparrow A \to hX$ including the gluon saturation effect with the recent preliminary experimental data from the PHENIX and STAR collaborations at RHIC. While we find overall reasonable agreement with the STAR data, our results indicate that the strong nuclear suppression of the asymmetry $A_N\sim A^{-1/3}$ observed by the PHENIX collaboration cannot be explained within the present understanding of this problem.


I. INTRODUCTION
Transverse single spin asymmetries (SSAs), as measured in collisions of an unpolarized probe with a transversely polarized proton, are traditionally a venue to understand the spin structure of the proton [1][2][3]. For inclusive hadron productions at high-P h⊥ , SSA is computed from perturbative QCD where it becomes a probe of collinear twist-3 distributions. Recent measurements at RHIC considered collisions of polarized protons on nuclear targets and so a completely new interplay between spin physics and the physics of gluon saturation becomes a reality [4][5][6][7][8][9][10][11]. This is especially so, as gluon saturation is important in the forward region of the produced hadron where SSA is the largest. Both the STAR [12] and the PHENIX collaborations [13] reported on preliminary results on SSA in p ↑ A → hX in addition to p ↑ p → hX. The PHENIX collaboration found a striking nuclear suppression A N ∝ A −1/3 with the mass number A in their preliminary data sets. On the other hand, the STAR collaboration did not find any significant nuclear effect in their data. The two data sets are not necessarily in contradiction to each other, as they are collected for different kinematics. However, the difference in kinematics is actually not very large, and both data are sensitive to the small-x region of the nucleus target. Therefore, it remains a challenge for theorists to explain both data consistently in a single framework.
At first sight, the suppression A N ∝ A −1/3 seems consistent with the prediction of k ⊥ -factorization approaches [4,5] which include the gluon saturation effects [14,15] in the target nucleus. However, the k ⊥ -factorization does not apply to this process, and a more proper treatment based on the collinear or 'hybrid' [9] factorization has identified two contributions with different scaling behaviors [16][17][18]. The recent fits of the p ↑ p data [19,20] indicate that the O(A 0 ) terms are dominant, so the PHENIX result is actually surprising. Furthermore, the suppression is observed at relatively high values of the hadron P h⊥ where one does not expect to see strong nuclear effects at RHIC energies. In view of this, it is premature to link the PHENIX finding with the gluon saturation effect.
In this paper we quantitatively address this problem by numerically computing the SSA in p ↑ p → hX and p ↑ A → hX using the formulas derived in [16,17]. We then compare our results with the preliminary STAR and PHENIX data. Following [19,20], we assume that the twist-3 fragmentation contribution is the main cause of SSA in this channel. As for the nucleus, we use the solution of the running coupling Balitsky-Kovchegov (rcBK) equation [21,22]. The main result, presented in Sec. III shows an overall satisfactory agreement with the preliminary STAR data. On the other hand, we were not able to confirm the nuclear suppression as seen in the PHENIX preliminary results. To the contrary, our results show no nuclear dependence for the PHENIX kinematics, even though we include the saturation effect of the target. We investigate the reason of this failure and discuss what extra contributions are needed to fix this problem.

II. CROSS SECTION FORMULAS
Our starting point are the formulas for the spin independent pA → hX cross section and the fragmentation contribution to the spin dependent p ↑ A → hX cross section within the CGC framework. The spin independent cross section [23] is given as where P h⊥ and y h are the hadron transverse momenta and rapidity, respectively, with x q,g = P h⊥ z √ s e ±y h and z min = P h⊥ √ s e y h and we sum over quark flavors as a . The function F (x, k ⊥ ) is the gluon dipole distribution defined as where Y = log(1/x), R A is the nuclear radius and U (x ⊥ ) is the fundamental Wilson line with . . . in the third line denoting the color average. In Eq.
(1) f a (x, Q 2 ) is the unpolarized parton distribution function and D h/a (z, Q 2 ) is the unpolarized hadron fragmentation function evaluated at the scale Q 2 = P 2 h⊥ . The spin dependent cross section comes from the quark-gluon-quark contribution, the twist-3 fragmentation contribution and the triple gluon contribution. In this work we consider only the fragmentation contribution as it is the dominant source of SSA in this channel [19,20]. We start from the main formula (see Eq. (46) in [17]) where M is the proton mass and S i ⊥ is the proton spin. Here h a 1 (x, Q 2 ) is the quark transversity distribution, while Imẽ h/a (z, Q 2 ) and ImÊ h/a F (z , z, Q 2 ) are the hadron twist-3 fragmentation functions. In the next step we approximate is a monotonically dropping function of k ⊥ . With this approximation it is possible to make use of the following relations for the twist-3 fragmentation functionŝ to eliminate the terms containing the z integral over ImÊ ) is yet another twist-3 fragmentation function. The notation used in this work relates to the notation in Ref. [20,24] as named as intrinsic, kinematical, and the dynamical twist-3 fragmentation functions. The equations (4) are known as the QCD equation of motion relation [25,26] and the Lorentz invariance relation [24], respectively. Using (4), Eq. (3) becomes In the following we will numerically compute the SSA defined as where in the numerator (denominator) we have the spin dependent (independent) cross section defined by Eq. (6) (Eq. (1)). We adopt the convention by which are azimuthal angles of the outgoing hadron (spin). When the incoming proton is pointing in the +z direction, and with its spin pointing in the y direction, ∆σ(↑) (∆σ(↓) = −∆σ(↑)) is the cross section for the hadron emission in the +x (−x) direction, or left (right) direction. This explains the "left -right" convention which is also used by STAR and PHENIX.
The nuclear effects are contained in the dipole function F (x, k ⊥ ) and especially the first term in the spin dependent cross section (6) depends on the derivative of the dipole. In the saturation regime (k ⊥ Q S ), where Q S is the saturation scale, we would typically get dF/dk ⊥ ∼ k ⊥ F/Q 2 S . Since the spin independent cross section (1) goes as ∼ F (x, P h⊥ /z), we find, for this particular term, A N ∼ Q −2 S , leading to A N ∼ A −1/3 for the nuclei. Although not immediately obvious, the second term of (6) also scales as A −1/3 [17]. From a quantitative perspective it is important that the saturation scale in the nuclei scales as (Q A S ) 2 = cA 1/3 (Q p S ) 2 (Q p S is the saturation scale in the proton) where an additional proportionality factor c < 1 [27] (in the numerical calculations we will use c = 0.5) will inhibit the overall magnitude of the nuclear suppression. On the other hand, when k ⊥ Q S we are in the perturbative regime where the dipole distribution has a characteristic dependence F ∼ Q 2 S /k 4 ⊥ and so dF/dP h⊥ ∼ Q 2 S /P 5 h⊥ . The same Q S -dependence is found also for the second term in (6): F/P h⊥ ∼ Q 2 S /P 5 h⊥ and so in the perturbative limit the nuclear dependence drops out in the ratio.

III. CALCULATION SETUP AND NUMERICAL RESULTS
In this Section we first explain all the details of our calculation and then we numerically compute SSA and compare with the available preliminary data from STAR and PHENIX. We will often be using the Feynman-x variable: For the dipole gluon distributions F (x, k ⊥ ) (2), we use the numerical solution of the running coupling Balitsky-Kovchegov (rcBK) equation [21,22] from [27]. We take the McLerran-Venugopalan (MV) initial condition at Y 0 = log(1/x 0 ), where x 0 = 0.01 as Here Q p,A S,0 is initial saturation scale parameter for the proton and the nuclei, γ is the anomalous dimension and Λ the IR cutoff. We use two different parameter sets for the initial condition. Labeling them as set MV and set MV γ , the model parameters are: MV: γ = 1, (Q p S,0 ) 2 = 0.2 GeV 2 and Λ = 0.241 GeV, MV γ : γ = 1.119, (Q p S,0 ) 2 = 0.168 GeV 2 and Λ = 0.241 GeV. For the nuclei, we use the relation (Q A S,0 ) 2 = cA 1/3 (Q p S,0 ) 2 with c = 0.5, as mentioned previously. For the twist-2 distribution functions f a (x, Q 2 ), we use the central CTEQ10 set [28]. For the twist-2 pion and kaon fragmentation functions D π/a (z, Q 2 ) and D K/a (z, Q 2 ) we use the central DSSV set [29]. The transversity distribution h a 1 (x, Q 2 ) and the twist-3 pion fragmentation functions Imẽ π + /a (z, Q 2 ) are obtained by solving their respective evolution equations numerically with the initial condition determined in [30]. The twist-3 kaon fragmentation function is obtained from [31]. In both cases we employ the Wilczek-Wandzura approximationê h/ā 1 (z, Q 2 ) = zImẽ h/a (z, Q 2 ) as in [30].
In Fig. 1 we show the numerical results of our computation for A N in p ↑ p → π 0 X and p ↑ Au → π 0 X as a function P π⊥ for several values of x F as compared to the preliminary STAR data [12]. We have used (Q Au S,0 ) 2 = 3(Q p S,0 ) 2 . The full (dashed) lines correspond to the calculation with the MV (MV γ ) model. The shaded band comes from the uncertainty in the transversity and the twist-3 fragmentation function from the analysis in [30]. In this case, we have calculated A N with the MV model. As a consistency check, we have also computed A N in p ↑ p using the collinear gluon PDF for the unpolarized proton (as was done in [20]). The result is in good agreement with the one from the rcBK solution shown in Fig. 1. 1 While the central results, given by the full and dashed lines on Fig. 1, seem to compare well with the preliminary STAR data within its experimental uncertainties, the shaded bands reflect a large theoretical uncertainty in the extraction of the transversity and the twist-3 fragmentation function. Nonetheless, there is a valuable point to be made here regarding the nuclear dependence of A N . The nuclear suppression of A N for p ↑ Au relative to A N in p ↑ p : AN for p ↑ p → π 0 X and p ↑ Au → π 0 X as a function of P π⊥ at √ s = 200 GeV versus preliminary STAR data [12]. The full (dashed) lines are a calculation using the MV (MV γ ) model where we used (Q Au S,0 ) 2 = 3(Q p S,0 ) 2 as well as the central values for the transversity and the twist-3 fragmentation functions, while the shaded band reflects the uncertainty in extraction of both of these quantities according to [30].
that is visible in the low P π⊥ region is not only the most prominent feature of our calculation, but also quite robust, being of a similar magnitude for the central results as well as for the shaded regions. Even though nuclear effects are not clear from the data itself, the overall magnitude of the nuclear suppression from both the central results and the shaded regions are still consistent with the data uncertainty.
In the preliminary PHENIX data set [13], covering a kinematics range x F ≤ 0.12, A N is measured in p ↑ p → h + X, p ↑ Al → h + X and p ↑ Au → h + X where h + is a mixture of outgoing π + and K + . The nuclear dependence of A N in the PHENIX data is most noticeable for the largest measured x F = 0.12 where also an average P h⊥ = 2.9 GeV is measured. At these kinematics, y h = 2.13 leading to x q > 0.12, x g > 0.0017 and so it is reasonable to apply the forward CGC formulas to compute A N . In Fig. 2 we show a comparison of our results for the nuclear dependence of A N in p ↑ A → h + X 2 to the preliminary PHENIX data, for the kinematics point x F = 0.12, P h⊥ = 2.9 GeV. In the PHENIX result A N clearly drops with the increase in the atomic number A, and this is consistent with the behavior A N ∼ A −1/3 . However, our current numerical results show virtually no A-dependence. The reason is clear: P h⊥ = 2.9 GeV is too hard to be sensitive to the saturation scale which is Q Au S ∼ 0.9 GeV for the PHENIX kinematics in the model used here.
To elaborate on this point, let us make an extreme assumption that only the first term of (6) is important. For P hT < Q A S , this term is expected to give the scaling A N ∼ A −1/3 and this is demonstrated in Fig. 3 where we plot the double ratio as a function of k ⊥ for several values of x using (Q Au S,0 ) 2 = 3(Q p S,0 ) 2 . Close to the initial condition the distribution is nearly Gaussian, and hence the ratio has a plateau in the low k ⊥ region at R N (Q p S ) 2 /(Q A S ) 2 1/3. For high k ⊥ , R N → 1 as a consequence of the perturbative tail. Going lower in x via the rcBK equation, the peak position  6)) as a contribution to AN . The shaded blue band takes into account the uncertainty of the transversity and the twist-3 fragmentation function according to [30]. moves toward the high-k ⊥ region and the plateau shrinks-already for x ∼ 0.001 the value of 1/3 is reached only for very small k ⊥ . For k ⊥ as large as 2.9 GeV, R N never deviates significantly from unity. Therefore, even in this extreme scenario A N does not have nuclear dependence. This is also illustrated by the thin curves in Fig. 2 where we computed A N including only the derivative term in the numerator. Note that we have only included the fragmentation contribution, but it is clear that adding the contribution from the twist-three quark-gluon distributions [16] will not help resolve this issue.

IV. CONCLUSIONS
We have made a numerical computation of SSA in p ↑ p and p ↑ A in the forward region including the gluon saturation effect of the nucleus. Using the current state-of-the-art twist-3 fragmentation functions and the dipole gluon distribution, we compared our results to the preliminary STAR and PHENIX data. While the saturation based description seems to describe well the STAR data, it fails to explain the scaling A N ∼ A −1/3 observed by the PHENIX collaboration.
According to the result of [16,17], a strong nuclear suppression of A N is possible only if the ∼ A −1/3 terms dominate over the ∼ A 0 terms, and if one looks at P h⊥ less than Q S . The recent fit of the pp data [20] suggests that the first condition does not hold 3 , and the second condition is also violated by the high value of P h⊥ measured by the PHENIX collaboration. This makes the PHENIX result all the more striking. It is also puzzling that there seems to be a sudden change in the behaviors of A N between x F = 0.2 (the lowest value measured by the STAR collaboration) and x F = 0.12 (highest value measured by the PHENIX collaboration).
This may call for alternative mechanisms of SSA around x F ∼ 0.1 whose nuclear dependence come from a different source. Indeed the region 0.1 x F 0.2 might be special-it is roughly the 'threshold' region where A N starts to grow. Thus the value of A N itself is very small in this region and a small effect can cause a large numerical impact. Perhaps one should include not only the q + g channel (as we did in this paper), but also the g + g channel. Since there is no 'gluon transversity' distribution, the fragmentation contribution is absent in this channel, but the collinear three-gluon, or 'odderon' contribution comes into play [7,10,[32][33][34]. A precise evaluation of SSA including these effects appears to be a challenging task.