Forward particle production in proton-nucleus collisions at next-to-leading order

We consider the next-to-leading order (NLO) calculation of single inclusive particle production at forward rapidities in proton-nucleus collisions and in the framework of the Color Glass Condensate (CGC). We focus on the quark channel and the corrections associated with the impact factor. In the first step of the evolution the kinematics of the emitted gluon is kept exactly (and not in the eikonal approximation), but such a treatment which includes NLO corrections is not explicitly separated from the high energy evolution. Thus, in this newly established"factorization scheme", there is no"rapidity subtraction". The latter suffers from fine tuning issues and eventually leads to an unphysical (negative) cross section. On the contrary, our reorganization of the perturbation theory leads by definition to a well-defined cross section and the numerical evaluation of the NLO correction is shown to have the correct size.


Introduction
When doing perturbation theory in Quantum Mechanics, e.g. when calculating the energy change in the states of a system due to a small perturbation, normally it suffices to sum a few terms. That is, E n = k 0 k=0 g 2k E (k) n , with g 2 the size of the perturbation, E (k) n fixed (as they depend only on a few parameters of the unperturbed Hamiltonian) and k 0 an integer chosen according to the desired accuracy. In an interacting Quantum Field Theory, like QCD, very often one is interested in a cross section, which perturbatively reads σ = σ 0 +g 2 σ 1 +g 4 σ 2 · · · . Now things can be different since the coefficients σ k can also depend on the momenta of the outgoing particles: although the QCD coupling α s = g 2 /4π is small at short distances, for some momenta one may have α s σ k+1 ∼ σ k . Then the fixed order series does not provide accurate information and one must resum to all orders in α s in such kinematical domains.
Even one of the most basic emissions in QCD, that of a gluon from a parent parton (quark or gluon), exhibits such features. Starting from the amplitude in Fig. 1.a, we find that the differential probability under consideration in the double logarithmic limit 1 reads with C R the Casimir for the representation of the parent parton. Say we want to measure the number of gluons with a very small longitudinal momentum fraction x. Then, at order n+1 in perturbation theory all intermediate gluons with strongly ordered longitudinal momenta as in Fig. 1.b are integrated over, leading to a large "coefficient" proportional to ln n (1/x). When x is sufficiently small, this coefficient becomes comparable to (or larger than) α n s , and twe must sum over all n, i.e. over all the diagrams with an arbitrary number of strongly ordered gluons. We find that the result exponentiates and, for example, the integrated gluon distribution (the number of gluons per unit rapidity and with transverse momentum less than Q 2 ) grows like Such a growth cannot continue to arbitrary small-x, since gluons start to overlap and the interaction among them becomes strong, even though the coupling is still weak at the transverse scale of interest. "Gluon saturation" will occur at a scale for which the number of gluons per unit phase space becomes of order 1/α s , that is when xG(x, Q 2 s )/Q 2 s πR 2 ∼ 1/α s , where R is the size of the hadron under study. Then the saturation momentum for a large nucleus of atomic number A reads with Q 2 0 being non-perturbative, while λ s can be determined perturbatively. For small-x and/or large A, the saturation scale is much larger than Λ QCD and indeed weak coupling techniques are applicable. Still, apart from the resummation of large logarithms, one eventually needs to deal with non-linear dynamics due to saturation. A QCD effective theory for describing such phenomena is the CGC [1].
The soft part of the hadronic wavefunction can be probed, for example, in proton-nucleus collisions in the forward region of the proton. A quark collinear with the proton and with a moderate to large momentum fraction x p = p + /Q + , multiply scatters with the soft gluons of the nucleus that carry a fraction X g = p − /P − , and picks up a transverse momentum k ⊥ so that it moves at an angle θ with the collision axis. Here, Q + is the longitudinal lightcone momentum of the proton which moves in the positive direction, while P − is that of the nucleus moving towards the negative one. A simple kinematics exercise shows that the quark rapidity in the COM reads η = − ln tan(θ/2), while x p = k ⊥ e η / √ s and X g = k ⊥ e −η / √ s, with s = 2Q + P − the COM energy squared. Then it is obvious that the forward region θ ≪ 1 maps to small-X g . We also note that in the infinite momentum frame of the nucleus, which shall be used later on, one has X g = k 2 ⊥ /x p s. Resumming both [ᾱ s ln(1/X g )] n andᾱ s [ᾱ s ln(1/X g )] n terms, withᾱ s = α s N c /π and N c the number of colors, in the presence of saturation, one has been able to give the cross section for single inclusive particle production at NLO in the CGC framework [2,3]. The problem is that such a cross section when evaluated numerically turns out to be negative (although it is positive at LO) for values k ⊥ ≳ Q s [4]. Here we shall review the process under consideration, and provide for the solution in the aforementioned problem [5].

Particle production at LO and the BK equation
In our analysis we shall omit the fragmentation of the outgoing forward quark into hadrons, since it does not play any role in the problem under discussion and as it is anyway straightforward to insert it at any time in the calculation. At LO, the collinear with the proton quark scatters multiply off the target field A µ , which can be strong for transverse momenta around (or smaller than) Q s . The quark moves having its transverse coordinate fixed at x, while its color undergoes a precession and such an eikonal interaction is given by the Wilson line where P stands for path ordering along the light-cone time x + of the quark and t a is in the fundamental representation. The direct amplitude (DA) M i j (k) is just the two-dimensional Fourier transform of V i j and then the single inclusive production for a particle (here the quark) with transverse momentum k and rapidity η is obtained by multiplying with the complex conjugate amplitude (CCA) (see Fig. 2.a) and by averaging/summing over colors. It reads where x p q(x p ) is the quark pdf in the proton, r = x − y and S (r; X g ) is the elastic S -matrix for the qq dipole (x, y) to scatter off the nucleus 2 (see Fig. 2.b). This S -matrix is bilinear in the Wilson lines (as we have taken the product of the DA with the CCA), and is given by It is evaluated at the target scale X g , which is determined from the kinematics, and resums the aforementioned powers ofᾱ s ln(1/X g ). In the problem at hand, we have already assumed in Eq. (5) that the nucleus is homogeneous so that S does not depend on the impact parameter b = (x + y)/2. When expanding the gauge field A µ to second order, i.e. when neglecting the multiple scattering, the cross section becomes proportional to the nuclear gluon pdf as it should. We do not elaborate on the details of the averaging in Eq. (6), but it suffices to say that it amounts to specifying the initial condition of an equation that we now discuss. The easiest way to implement the resummation of the large logarithms is to write the BK equation [7,8] which describes the evolution of the elastic dipole S -matrix. For our purposes,  it is more convenient to present such an equation in its integral form which reads where X(x) = X g /x is the longitudinal momentum fraction in the target in the arbitrary intermediate evolution step. The kernel in Eq. (7) is called the dipole kernel [9] as it stands for the transverse dependence of the probability for the dipole r to split into z and r − z. S (r, X 0 ) is the tree level contribution, the one in the absence of any QCD evolution, and typically it is given by the McLerran-Venugopalan (MV) model [10] or some variant of it. Such an equation is better understood in a frame where almost all of the evolution is associated with the nucleus, apart from the emission of the primary gluon (the one closest to the dipole r), cf. Fig. 3. It is clear that the solution to Eq. (7) unitarizes, i.e. S → 0 when X g → 0 for a fixed dipole size r ⊥ = |r| (as we assume that S depends only on the magnitude of r). In the current setup, it is natural to define the saturation momentum as the line in the r ⊥ -X g plane along which the S -matrix is fixed to a value of order (but of course smaller than) 1, e.g.
Indeed one finds that Q 2 s (X g ) grows as an inverse power of X g , although the growth is too strong. Using NLO BK dynamics and resumming collinear logarithms which can get large, the growth is tamed [11][12][13] and comes to a good agreement with the phenomenology. Moreover, eventually one is able to reasonably describe the relevant data in d-Au collisions using Eq. (7) with running coupling and inserting K-factors to account for the remaining higher order corrections [14,15].

Particle production at NLO
As depicted in Fig. 4.a, NLO corrections to the single inclusive particle production, can be separated in two parts. First, one may have two adjacent gluons in the "ladder" which have fractions such that x n ∼ x n+1 ≪ 1, or one may have a loop. These are corrections associated with the evolution of the color dipole and can be taken into account by using the NLO BK equation. We shall not deal with them, since in principle they can be implemented after collinear resummations have been done without creating any additional issues. Second, one has a NLO correction when the primary gluon is not emitted eikonally, i.e. when its fraction x is of order (but smaller than) 1. Thus, by treating exactly the kinematics one can include this  impact factor correction (still, one must notice that the scattering with the nucleus remains eikonal and is given by Wilson lines) into the calculation of the cross section. Then, by extending the LO result in Eqs. (5) and (7), it is not hard to see that at NLO we have the sum of the two terms in Fig. 4.b. In fact there are two NLO contributions, one proportional to N c and one to C F , and we focus on the former which reads The real J and virtual J v pieces (implicitly proportional toᾱ s ) are known [2,3,16], S stands for the Fourier transform of S , and 1 − ξ = x is the fraction of the primary gluon. In a compact notation, where we omit the quark pdf, suppress the transverse coordinates and trivial prefactors, and portray collectively the scattering of the qqg system, we can write Eq. (10) is our main result and being the sum of the tree level contribution and a positive term should not pose any problems if no further approximations are done [5]. Indeed, this will be manifest in the numerical results to be shown in the next section.

k ⊥ -factorization: fine-tuning and negativity issues
Eq. (10) is not in the traditional form of k ⊥ -factorization. This is so, because the kernel K depends not only on the transverse coordinates, but also on x, and moreover, the S -matrices of the qqg system are evaluated at the "floating" scale X(x). Now the task is to see how we can arrive at a k ⊥ -factorized expression. First notice that the LO result is recovered, as it should, by letting in Eq. (10) K(x) → K(0): indeed, the two terms in the r.h.s. combine to give S (k, X g ). Then, we are free to add and subtract this LO piece to the r.h.s. of Eq. (10) and by doing so we get the "subtracted" expression This is a mathematically correct step, however it is a dangerous one, since there is reshuffling of a very large contribution between the two terms. In order to obtain the correct result from Eq. (11), one must be extremely precise on the numerical evaluation of the two terms [5]. Two more steps, which now involve approximations, are required to bring Eq. (11) in a k ⊥ -factorized form. Since K(x) − K(0) clearly vanishes as x → 0, the integrand on the r.h.s. is not anymore dominated by small values of x, and therefore at NLO accuracy we can let S qqg k, X(x) → S qqg k, X g and furthermore let X g → 0 in the lower limit of the integration. Then we arrive at [5] dN k ⊥−fact dηd 2 k = S(k, X g ) +ᾱ s which is what "standard" k ⊥ -factorization would give: the S -matrices are evaluated locally at X g , while the kernel depends only in the transverse coordinates (since, in principle, we can now integrate over x). We repeat that, while the "unsubtracted" and "subtracted" expressions (10) and (11) are mathematically identical, Eq. (12) is not equivalent with them. Due to the approximations performed, it mistreats the kinematics at small-x and leads to a negative cross section. Indeed, all this is confirmed in the numerical evaluation of all the aforementioned formulas for the cross section as shown in Fig. 5. We also note that, even though Eq. (10) and Eq. (11) agree very well with each other, the subtracted expression exhibits some fluctuations for high k ⊥ as a result of the large reshuffling. One sees that the NLO correction is negative and with a large but under control magnitude: around 50% at k ⊥ ∼ 5 GeV [17].
The discussion so far has been restricted to using a QCD coupling fixed to a moderately small value. It is not trivial at all to extend Eq. (10) with running coupling included. The dependence onᾱ s there is in two places: the implicit BK-dependence in S qqg k, X(x) and the explicit impact factor one in front of the integration. Regarding the former, that is the coupling along the small-x ladder, one first solves the BK-equation in coordinate space with the smallest dipole prescription [18], i.e. usingᾱ s (r min ) with r min the size of the smallest dipole among the three involved in the splitting process in Eq. (7) and then one transforms the results to momentum space. Regarding the latter, that is the coupling associated with the primary gluon, it is natural to useᾱ s (k ⊥ ) [5], since k ⊥ is the momentum of the "measured" outgoing quark, and all this is leading again to well-defined results. In order to have a more homogeneous description, one may try to move the coupling inside the integration over the transverse coordinates [17]. Without going into much detail here, we find that when the primary gluon has a moderate x and is soft in the transverse space, the scale in the coupling must be set by the "daughter dipole" [19], which is the largest one in the kinematics under consideration. A smallest dipole choice here would lead to heavily unstable results [17], since it generates a fake potential leading to unphysical contributions to the scattering [19]. Moreover, for the NLO terms in the C F -sector a coordinate space prescription for the coupling either cannot be implemented, or leaves uncanceled some spurious longitudinal logarithms [19], so eventuallyᾱ s (k ⊥ ) is the only viable choice.

Conclusion
Single inclusive particle production at forward rapidities in proton-nucleus collisions is a process suited to study gluon saturation. Early results at NLO in the CGC effective theory led to a negative cross section and here we have shown that this was the result of approximations which were necessary to arrive at an expression consistent with k ⊥ -factorization. Thus, we have developed a generalized factorization procedure (which involves a longitudinally dependent kernel and is non local in rapidity) that gives by definition a well-defined cross section [5]. As expected, the NLO correction is negative and reduces the LO result by ∼ 50% when the outgoing quark (that fragments to a hadron) has a transverse momentum of a few GeV.