On the use of a running coupling in the NLO calculation of forward hadron production

We address and solve a puzzle raised by a recent calculation [1] of the cross-section for particle production in proton-nucleus collisions to next-to-leading order: the numerical results show an un- reasonably large dependence upon the choice of a prescription for the QCD running coupling, which spoils the predictive power of the calculation. Specifically, the results obtained with a prescription formulated in the transverse coordinate space differ by one to two orders of magnitude from those obtained with a prescription in momentum space. We show that this discrepancy is an artefact of the interplay between the asymptotic freedom of QCD and the Fourier transform from coordinate space to momentum space. When used in coordinate space, the running coupling can act as a fictitious potential which mimics hard scattering and thus introduces a spurious contribution to the cross-section. We identify a new coordinate-space prescription which avoids this problem and leads to results consistent with those obtained with the momentum-space prescription.


I. INTRODUCTION
Particle production in proton-proton or proton-nucleus collisions at forward rapidities and semihard transverse momenta of a few GeV represents an important source of information about the small-x part of the nuclear wave function, where gluon occupation numbers are high and nonlinear effects like gluon saturation and multiple scattering are expected to be important. To optimize our extraction of this physical information from the experimental data at RHIC and the LHC, reliable theoretical predictions are of the utmost importance. For such semihard momenta, the QCD coupling is only moderately weak; hence, the perturbative calculations of the relevant cross sections must be pushed up to at least next-to-leading order (NLO) accuracy.
In the recent years, the color glass condensate (CGC) effective theory [2] (the natural framework for such calculations within perturbative QCD) has indeed been promoted to NLO accuracy. This includes both the equations describing the high-energy evolution of the scattering amplitude in the presence of nonlinear phenomena-the B-JIMWLK. 1 hierarchy of equations [3][4][5][6][7][8][9] and its mean field approximation known as the Balitsky-Kovchegov (BK) equation [3,10]-and the impact factors describing the coupling between a dilute projectile and the dense gluon distribution from the nuclear target. However, the first NLO calculations met with serious difficulties (instability of the NLO BK equation, negative cross section for particle production, huge scheme dependence in the choice of a prescription for the running of the coupling) that were progressively understood and overcome.
Specifically, the NLO version of the BK equation [11][12][13] turned out to be unstable [14,15], due to the presence of large and negative NLO corrections enhanced by double collinear logarithms. Two independent solutions have been proposed to solve this problem-the enforcement of a kinematical constraint 2 [17] and a collinearly improved version of the BK equation [18,19]-which refer to two different methods for resumming the double collinear logarithms to all orders. After restoring the full NLO accuracy (by adding the remaining NLO corrections, in particular, those expressing the running of the coupling) on top of the collinear improvement [18,19], the numerical solutions to the BK equation were found to be stable and physically meaningful [20].
Furthermore, the first NLO calculations of single-inclusive particle production in proton-nucleus collisions at forward rapidities [21][22][23][24][25][26][27] led to a cross section which suddenly drops and becomes negative when increasing the transverse momentum of the produced hadron (but still within the semihard region where the formalism is supposed to apply). These NLO calculations used the so-called "hybrid factorization" [28] together with the NLO correction to the impact factor originally computed by Chirilli, Xiao, and Yuan [21,22] (see also [29] for an alternative calculation). The word "hybrid" refers to the fact that the projectile proton and the nuclear target are treated on a different footing. The proton, which is dilute, is treated within collinear factorization, that is, by using the standard, "integrated," parton distributions which also appear in the study of deep inelastic scattering. The dense nucleus, on the other hand, is treated within the spirit of k T factorization, that is, in terms of "unintegrated" gluon distributions which describe the distribution of the (soft) gluons in both longitudinal and transverse momenta, including its nonlinear evolution with increasing energy. The "impact factor" encodes the coupling between a collinear parton from the incoming proton and the soft gluons in the nucleus. Within k T factorization, it is assumed that this impact factor can be computed at a fixed rapidity separation between the projectile and the target, and hence, it can be separated from the high-energy evolution. This is however just an approximation: starting with NLO, the calculation includes gluon fluctuations which span the whole rapidity phase space and thus probe the target gluon distribution within an extended range in rapidity.
As explained in Ref. [30], the problem with the negativity of the NLO cross section mainly comes from enforcing such a local (in rapidity) separation between the evolution and the NLO impact factor (the "rapidity subtraction"). This difficulty is therefore expected to be generic. And indeed, as recently shown in Ref. [31], a similar problem appears in the context of deep inelastic scattering (DIS), when using the "dipole factorization" (a version of k T factorization appropriate for DIS at high energy) together with the NLO impact factor from Refs. [32,33].
An additional source of difficulty, that will be our main focus in this paper, is the mismatch between the running coupling prescriptions used in the calculation of the NLO impact factor and in the high-energy evolution (the BK equation), respectively. On one hand, the BK equation is most naturally formulated and solved in the transverse coordinate space. Indeed, the coordinate representation allows for simple implementations of the eikonal approximation and of the unitarity constraint on the dipole scattering amplitude. On the other hand, the cross section, hence the NLO impact factor, must be ultimately expressed in terms of the transverse momentum of the produced particle, which is a measurable quantity. 3 This makes it natural to use a coordinate-space scale (a dipole size) for the running coupling in the solution to the BK equation, but a momentum-space scale (a transverse momentum) in the one-loop calculation of the NLO impact factor. (The precise prescriptions that we shall use in practice will be specified later on.) Albeit not a fundamental problem in itself-the use of various running coupling (RC) prescriptions in different sectors of the calculation is consistent with the overall NLO accuracy and should be viewed as a part of the scheme dependence inherent in the formalism-this mismatch amplifies the negativity problem of the NLO cross section [30]. In fact, it may introduce such a problem by itself, even when the separation between the leading-order result, and the NLO corrections is made by properly keeping the nonlocality of the one-loop calculation in rapidity.
To avoid such problems altogether, Ref. [30] proposed an "unsubtracted" factorization scheme which avoids the subtraction of the leading-order evolution from the full NLO result. This factorization preserves the actual skeleton structure of the original Feynman graphs: the "impact factor" is now given by a full one-loop graph in which the soft (smallx) part of the gluon fluctuation contributes to the (leadingorder) BK evolution, whereas its hard [x ∼ Oð1Þ] part contributes to NLO. This "unsubtracted" formulation yields a positive-definite cross section by construction, including when using different RC prescriptions in the impact factor and in the high-energy evolution, as above explained. This has been explicitly verified by the numerical calculations in [1], which also found that the NLO corrections are negative and significantly large-they reduce the cross section by up to 50% with respect to the respective prediction at the leading order (LO). A similar behavior has been recently observed in the NLO calculation of the DIS structure functions [31]: once again, the "unsubtracted" version of the factorization leads to positive results.
Reference [1] has also considered an alternative formulation of the NLO calculation, which differs from the "canonical" one in the presence of RC corrections, in order 2 See also Ref. [16] for a generalization of the kinematical constraint to the Langevin formulation of the JIMWLK equation. 3 The transverse momentum representation of the impact factor is also convenient in view of the subtraction of the collinear divergences from the one-loop calculation [22]; see the discussion in Sec. IV below.
to study the scheme dependence of the results. In this new formulation, the one-loop calculation of the impact factor is fully rewritten in coordinate space and the transverse momentum of the produced parton is introduced via a final Fourier transform. In this context, it is possible and also natural to use a coordinate-space argument for the RC also in the impact factor and thus have a unified treatment for the RC throughout the calculation. Clearly, the final results for the cross section need not be exactly the same in this scheme and the "canonical" one (which uses a momentum-space RC prescription in the impact factor), yet, they were expected to lie close to each other, because the Fourier transform roughly identifies momenta with the inverse of coordinates. So, it came as a real surprise (and also as a bad news for the reliability of the NLO formalism as a whole) when it turned out that the numerical results are dramatically different within these two schemes: the NLO corrections obtained in the new scheme have the opposite sign as compared to those in the "canonical" scheme, and they are tremendously larger-by 1 to 2 orders of magnitude, depending upon the final transverse momentum [1].
It is our main purpose in this paper to understand the origin of this puzzle and thus restore the predictive power of the NLO formalism. In our analysis, we shall trace this problem to the fact that the choice of a scale for the running of the coupling does not commute with the Fourier transform and show that this mathematical property is enough to explain the dramatic consequences observed in practice. This may look surprising in view of the fact that the scale dependence of the running coupling is quite weak-merely logarithmic. Yet, as we shall show, it is precisely this logarithmic decrease of the RC with increasing transverse momenta (or decreasing transverse sizes)-which is of course the expression of asymptotic freedom-which is responsible for the problem identified in our NLO calculation. Namely, when used in coordinate space, the size-dependent RC may effectively act as a fictitious "potential", which allows for hard scattering (due its logarithmic singularity in the limit of a zero size) and hence can transfer a very high momentum to the produced parton even in the absence of any physical scattering. This fictitious scattering produces a spurious component in the particle production, which at high transverse momenta decreases slower than the physical tail and hence dominates over the latter.
The physical interpretation of this technical argument becomes more transparent if one starts in coordinate space. In order to produce a very hard parton, with a transverse momentum much larger than the saturation momentum of the target, the produced particle must be accompanied by an unmeasured, recoil gluon, which carries the same transverse momentum but with the opposite sign. (It is precisely this recoil, or "primary" gluon, which is responsible for the oneloop contribution to the impact factor.) Soft primary emissions cannot contribute to the cross section, simply by momentum conservation. But when the one-loop calculation is performed in coordinate space, this seemingly trivial kinematical constraint is built in a rather nontrivial way. The loop correction is then dominated by gluons which are emitted far away from their parent parton and which physically correspond to soft emissions. Yet, such soft gluons do not contribute to the production of that parent parton-as expected from momentum conservation-because they are suppressed by the final Fourier transform. But this suppression can be spoilt by the scale dependence of the RC: that is, even soft gluons can formally contribute to the production of a hard parton because the necessary transverse momentum can be (incorrectly) provided by the RC. This unphysical contribution is the source of the surprising results obtained with a coordinate-space prescription in Ref. [1].
This discussion, developed in great detail in Sec. III, shows that one must be very careful whenever the use of a RC is accompanied by a Fourier transform. There are other situations where the interplay between the running of the coupling and the Fourier transform has been reported to lead to difficulties-e.g., to positivity violations for the "unintegrated gluon distribution" [34], defined as a Fourier transform of the dipole scattering amplitude. This last problem will however not occur in our subsequent calculations. Presently, it is not clear to us whether these two problems are somehow related at a deeper level.
This being said, we shall nevertheless be able to identify at least one RC prescription which is formulated in coordinate space and which circumvents the aforementioned problem of the "fake potential": this is the case where the argument of the RC is the transverse separation between the daughter gluons and their parent parton. This particular prescription quasicommutes with the Fourier transform, and in any case, it cannot act as an unphysical source of transverse momentum. Indeed, we shall see that NLO calculations using this prescription gives results which are extremely close to those of the momentum-space RC prescription, thus confirming that the overall scheme dependence of our NLO factorization is reasonably small.
The daughter-gluon prescription 4 looks particularly appealing in that it can be simultaneously used within the one-loop calculation of the impact factor and the BK equation, thus providing a coherent scheme for the ensemble of the NLO calculation. 5 Results obtained with such a unified prescription will be presented as well and found to be rather close to those predicted by the "canonical" mixed-space prescription (momentum-space RC for the impact factor, but 4 In the remaining part of this paper, we shall generally use the dipole picture of the BK evolution, as valid in the limit of a large number of colors; correspondingly, we shall refer to this prescription as the "daughter-dipole prescription". 5 Note however a limitation of this scheme in practice, to be discussed in more detail in Sec. IV: it cannot be used within the subclass of NLO corrections ("the C F terms") that are concerned with the subtraction of collinear divergences. Indeed, this last operation has so far been performed only in momentum space [22]. coordinate-space RC for the BK equation). However, one should keep in mind that the daughter-gluon prescription is not very well motivated in the context of the high-energy evolution-especially for the asymmetric problems where there is a large disparity between the characteristic transverse sizes of the target and of the projectile (see the discussion in Appendix A below). This is in particular the case for the forward particle production in proton-nucleus collisions, especially when the produced parton is relatively hard. In such a case, we would still recommend the use of a "mixedspace" running-coupling prescription.
The paper is organized as follows. Section II presents the NLO result for the quark multiplicity, as it will be used in this paper. Besides introducing the notations and the relevant formulas, this presentation will also give us an opportunity to anticipate some of the subtleties with the use of a running coupling. Section III is the main section of this paper. This is where we explain the origin of the puzzle with the coordinate-space RC prescription, as identified in Ref. [1]. We furthermore present an alternative RC prescription, still in coordinate space, which avoids this problem and can be meaningfully compared with the momentum-space prescription. The discussion is illustrated with numerical results. Section IV considers a special set of NLO corrections ("the C F terms") which require a different treatment, since for them the transverse coordinate of the primary gluon has already been integrated out. (More technical aspects of this discussion are deferred to Appendix B.) Finally, Sec. V summarizes our results and conclusions. For more clarity, we present in Appendix A the usual RC prescriptions used in relation with the BK equation together with a physical discussion.

II. SINGLE INCLUSIVE QUARK PRODUCTION AT NEXT-TO-LEADING ORDER
In this section, we shall summarize the next-to-leading order results for the quark multiplicity in proton-nucleus collisions at forward rapidities together with the underlying physical picture. Our main focus will be on the "unsubtracted" formulation of the hybrid factorization to NLO [30], but we shall also describe the original formulation by Chirilli et al. [21,22] and the associated problems. Finally, we shall explain some subtleties which occur when trying to consistently include the running coupling effects in the various sectors of the calculation [1,30].

A. The hybrid factorization to NLO
In the CGC effective theory [2], projectile partons are energetic enough so that their scattering off a large nucleus can be computed in the eikonal approximation: the only effect of the scattering is a color rotation described by a Wilson line in the appropriate representation of the color group SUðN c Þ. For example, for a right-moving quark, the Wilson line reads where g is the QCD coupling, x þ the light cone time of the projectile quark and x its transverse coordinate, A − a the relevant component of the color field of the nucleus, and t a the generators of SU(N c ) in the fundamental representation.
We shall be interested in forward particle production in proton-nucleus collisions and in particular, in the quark multiplicity, i.e., in produced quarks with transverse momentum k and rapidity η in the center of mass frame (COM). (Obviously, the respective gluon multiplicity will also contribute to the production of forward particles, but we shall not discuss its computation in the current work.) We choose the proton to be a right mover, hence the nucleus to be left mover, and work in the light cone gauge A þ a ¼ 0 in which the partonic structure of the proton (say, as encoded in the respective light cone wave function) is manifest. To compute the cross section for quark production, we use the so-called "hybrid factorization" generally assumed to apply to such a "dilute-dense" scattering-its validity has been so far demonstrated to next-to-leading order (NLO) in perturbation theory. In this factorization, the only information that we need to know about the dilute projectile (the proton) are the standard ("collinear") parton distributions-here, the quark distribution x p qðx p ; Q 2 Þ. In applications to phenomenology, one should also consider the fragmentation of the produced quark into hadrons, but here we shall ignore this final-state evolution, since it is irrelevant for the problems we would like to address.
Under these assumptions, the production of a quark at forward rapidity (i.e., for η positive and large) and to leading order (LO) in pQCD is the result of the (multiple) scattering between a "valence" (in the sense of large x p ) quark from the incoming proton and the dense gluon distribution of the nucleus. The quark, originally collinear with the proton and hence with no initial transverse momentum, acquires its final transverse momentum k through this scattering. The nuclear gluon distribution is dense since the typical gluons probed by this process have a very small value of the "minus" longitudinal momentum fraction X ≪ 1 (see below) and relatively small transverse momenta, q ⊥ ∼ Q s ðXÞ, hence their occupation numbers are large. Here, Q s ðXÞ is the nuclear saturation momentum at the longitudinal scale X relevant for the scattering.
To the same accuracy, the scattering between the quark and the color fields representing the gluons in the nuclear target can be computed in the eikonal approximation, cf. Eq. (2.1). Then the LO quark multiplicity is computed as 6 6 The resolution scale (Q 2 ) dependence of the quark distribution will be omitted throughout this paper; in practice [and for the relatively hard transverse momenta of interest here: where Sðk; XÞ is the "unintegrated gluon distribution" in the target, operationally defined as the Fourier transform of the S matrix describing the elastic scattering between a small color dipole and the nucleus, where in coordinate space we have The "color dipole" is made with a quark at x and an antiquark at y, which physically represent the produced quark in the direct amplitude and in the complex conjugate amplitude, respectively. In Eq. (2.3), r ¼ x − y is the dipole size, and we also assume, for simplicity, that the nucleus is quasihomogeneous in the transverse plane, so the S matrix is independent of the impact parameter b ¼ ðx þ yÞ=2. The average in Eq. (2.4) is taken over the color field A − a and with the CGC weight function which describes the high-energy evolution of our observable. Finally, the longitudinal momentum fractions for the quark, x p , and for the gluons, X g , which enter the formula (2.2) for the LO multiplicity are determined by the kinematics of the collision and of the produced quark, as follows: where k ⊥ ¼ jkj and ffiffi ffi s p is the COM energy. As anticipated, the forward kinematics (η positive and large) corresponds to the regime X g ≪ x p < 1.
Still at LO, one needs to resum to all orders the radiative corrections enhanced by powers ofᾱ s Y g , withᾱ s ≡ α s N c =π and Y ≡ lnð1=X g Þ, as generated by the high-energy evolution, i.e., by the successive emissions of soft gluons, which are strongly ordered in longitudinal momenta and whose effects should be computed in the eikonal approximation. To LO at least, this evolution can be associated with either the nucleus, or the incoming quark (or dipole), and in what follows, we shall mostly adopt the first point of view (although the viewpoint of dipole evolution will be useful too for some arguments 7 ). Sðx; y; XÞ is then the elastic amplitude for a "bare" dipole ðx; yÞ which scatters off the nuclear gluon distribution evolved in X from some initial value X 0 ∼ 1 [withᾱ s lnð1=X 0 Þ ≪ 1] down to the scale X g ≪ 1 of interest. To LO in pQCD and in the limit of a large number of colors, N c → ∞; this evolution is described by the Balitsky-Kovchegov (BK) equation [3,10] with z denoting the transverse coordinate of the soft gluon emitted in one step of the evolution. The rhs features the "dipole" version of the LO BFKL kernel. The initial condition Sðx; y; X 0 Þ for this equation should be computed with a suitable model for the nuclear gluon distribution at X 0 , like the McLerran-Venugopalan (MV) model [35,36].
At NLO, one needs to also include the "pure-α s " corrections, i.e., the radiative corrections of Oðα s Þ which are not enhanced by Y g . These can be divided into two classes: (i) NLO corrections to the high-energy evolution, i.e., to the kernel (more generally, to the structure) of the BK equation, and (ii) NLO corrections to the "impact factor", i.e., to the "hard" matrix element which describes the quark-nucleus scattering in the absence of any evolution, meaning already for X g ∼ X 0 . The LO impact factor describes the scattering between a bare quark collinear with the proton and the nucleus. At NLO, the wave function of the incoming quark may also contain a gluon with a "plus" longitudinal momentum fraction x (with respect to the parent quark). Hence, the NLO impact factor must describe the scattering between this quark-gluon system and the nuclear target. So long as x ≪ 1, this gluon emission can be computed in the eikonal approximation, and its effect can then be included in the LO evolution. So, strictly speaking, the NLO correction to the impact factor should only include relatively hard gluon emissions with x ∼ Oð1Þ, which must now be computed exactly (i.e., beyond the eikonal approximation). In practice though, it turns out that separating the LO evolution from the NLO correction to the impact factor is quite subtle (especially in the presence of running coupling corrections), as we shall see. For this reason, we shall prevent ourselves from doing such a separation, and mostly consider an "unsubtracted" expression for the NLO quark multiplicity [30] in which the two types of corrections are still mixed with each other.
As mentioned in the Introduction, our main goal in this analysis is to clarify the use of the running coupling in the calculation of the NLO multiplicity, notably in relation with the interplay between the high-energy evolution and the impact factor. For that purpose, we can ignore all the NLO corrections to the BK equation, except for those expressing the (one-loop) running of the coupling. That is, throughout this paper, we shall use the running-coupling version of Eq. (2.6), to be referred to as "rcBK", with various prescriptions for the argument ofᾱ s , to be later specified.
The NLO expression for the quark multiplicity can be decomposed into two pieces, one proportional to the gluon Casimir N c and the other one to the quark Casimir C F (see, e.g., [27,30] for details). In the formulation of Ref. [30], the first piece, proportional to N c , is the one which encodes both the LO evolution and the part of the NLO impact factor which is most relevant for the present discussion. So, over most of the subsequent developments we shall focus on this N c piece alone. (The C F -piece will be separately discussed, in Sec. IV.) After including the NLO corrections proportional to N c , the quark multiplicity reads [30] dN LOþN c The first piece in the rhs is the would-be "tree-level" result, proportional to the LO impact factor Sðk; X 0 Þ. The second term, expressed as an integral over ξ, encodes the dynamics of the (generally) noneikonal splitting of the incoming quark into a quark and the "primary" gluon, and the subsequent eikonal scattering of this quark-gluon system off the evolved nucleus. The two terms under the integral, J and J v , refer to real and, respectively, virtual gluon emissions, where we recall that a parton is considered as "real" if it appears as an on shell particle in the final state.
(Some illustrative Feynman diagrams will be shown in Sec. II C, in coordinate space.) The integration variable ξ represents the splitting fraction of the quark, that is, ξ ¼ 1 − x, with x the longitudinal momentum fraction of the primary gluon. Since the "plus" component of the scattering partons is not modified in the eikonal approximation, also the observed quark carries a fraction ξ of the parent quark. The latter has therefore a longitudinal fraction x p =ξ with respect to the projectile proton. On the other hand, in the "virtual" term, which is needed for probability conservation, the final quark has the same longitudinal momentum fraction, equal to x p , as the incoming quark. The fact that, when ξ ≠ 1, the J and J v terms in Eq. (2.7) are weighted by different values of the quark distribution function has therefore a trivial kinematical origin, but it will have important consequences in practice, as we shall see. Specifically, these two terms are given by [21,22]  In writing these expressions, we have used the large-N c limit in order to factorize the dipole S matrices appearing in the rhs.
On the other hand, the color algebra associated with the emission of the primary gluon and with the subsequent scattering of the quark-gluon system has been computed exactly, for a generic value N c . Hence, the overall factor N c (as encoded in the rescaled coupling constantᾱ s ¼ α s N c =π) can be trusted even beyond the large-N c limit. A similar discussion applies to the C F terms to be presented in Sec. IV.
Notice that J and J v have two sources of dependence upon ξ: one explicit in the transverse momentum kernels in Eqs. (2.8)-(2.9), which reflects the noneikonal nature of the splitting, and one implicit in the rapidity arguments XðξÞ of the dipole S matrices, which expresses the evolution of the nuclear gluon distribution from the initial scale X 0 down to the longitudinal scale XðξÞ probed by the effective projectile made with the quark and the primary gluon. This value XðξÞ depends upon ξ (and, in general, is larger than X g ) because the emission of the primary gluon "consumes" a rapidity interval Δy ¼ lnð1=xÞ, thus reducing the rapidity interval available to the evolution of the gluon distribution in the target, which now only goes up to Y g − Δy ¼ lnðx=X g Þ. These considerations suggest that where the last equality follows after also using Eq. (2.5). This turns out to be approximately correct, but only when the transverse momentum of the produced quark is large enough, k ⊥ ≳ Q s ðXÞ, with Q s ðXÞ the target saturation momentum (see [30] for details). This is indeed the most interesting kinematics for our purposes and the only to be considered in what follows. Notice that the constraints X g ≤ XðξÞ ≤ X 0 -as coming from energy-momentum conservation together with our choice for the initial condition for the BK evolution-imply indeed 0 ≤ ξ ≤ 1 − X g =X 0 , in agreement with the integration limits visible in Eq. (2.7). As anticipated, the quantities J and J v include both the LO evolution (in the sense of rcBK) and the NLO corrections to the impact factor, albeit in a rather subtle way. The part of the LO evolution associated with the emission of a soft primary gluon is generated by the region ξ ≃ 1 (i.e., x ≪ 1) of the integrals over ξ in Eqs. (2.8) and (2.9); in this region, one can set ξ ¼ 1 within the emission kernels, which is tantamount to working in the eikonal approximation. As for the subsequent, soft, gluon emissions responsible for the high-energy evolution over the remaining rapidity interval ln½X 0 =XðξÞ, they are of course encoded in the various S matrices, which are obtained by solving rcBK with initial condition Sðr; X 0 Þ and then making a Fourier transform to momentum space. Concerning the NLO corrections to the impact factor, these have two origins: the noneikonal (ξ ≠ 1 in the emission kernels) part of the primary gluon emission and the running coupling corrections to theᾱ s prefactor visible in Eqs. (2.8)-(2.9), which so far has been formally treated as a fixed coupling. The separation of the LO evolution from the NLO impact factor together with various prescriptions for the running couplings (in the primary vertex and the BK equation) will be discussed in the next sections.

B. Unsubtracted, subtracted, and CXY expressions
In this section, we shall decompose the NLO result in Eq. (2.7) into a leading order piece plus NLO corrections to the impact factor. Then we shall perform additional approximations in order to recover the original expression for the NLO quark multiplicity, by Chirilli et al. (CXY) [21,22], which features the "plus" prescription in the integral over ξ. In turn, this "plus" prescription is a variant of the "k T factorization", here applied to the unintegrated gluon distribution in the nucleus. Our main message from this discussion is that the "k T factorization", which is local in rapidity, is not equivalent to our general formula (2.7), but rather involves additional approximations which can be troublesome in practice. The subsequent manipulations are a priori valid at fixed coupling. Their extension to a running coupling is quite subtle and can bring additional complications, as we shall see in the next sections.
We start by introducing more compact notations, to be used only in this section. Namely, we rewrite Eq. (2.7) as where the definition of K is easily understood by comparing to Eq. (2.7). The first term in the rhs is the tree-level contribution or equivalently the initial condition, while the superscript in the second term stands for "unsubtracted" and will be shortly explained.
The LO contribution to the quark multiplicity should be recovered from Eq. (2.7) in the limit where the primary gluon emission is evaluated in the eikonal approximation. As already explained, this corresponds to replacing ξ → 1 within the kernels in Eqs. (2.8)-(2.9), as well as in the quark distribution accompanying the "real" emission, but not in the rapidity arguments XðξÞ of the various S matrices. This yields To see that this reproduces indeed the expected LO result from Eq. (2.2), we can use the identity which is (in compact notations) the Fourier transform of the integral version of the LO BK equation (2.6).
To separate leading from next-to-leading order contributions in Eq. (2.7), we subtract the LO result in the form of Eq. (2.12) and then add it back in its original form in Eq. (2.2). One thus finds where we have used the compact notation introduced in Eq. (2.11). Notice that the integrand in the "subtracted" piece has no singularity as ξ → 1 (i.e., x → 0), meaning that the integral over ξ develops no "small-X g " logarithm. This is in agreement with the fact that the longitudinal logarithm generated by a soft primary emission has been included in the evolution of the LO dipole Sðk; X g Þ.
Because of this subtraction, the integral in Eq. (2.14) is a pure-α s effect-a NLO correction to the impact factor associated with a relatively hard primary emission.
Clearly, Eqs. (2.11) and (2.14) are equivalent to each other, since related by exact manipulations. Notice however that in going from Eq. (2.11) to Eq. (2.14), one has added and subtracted a large term (the LO contribution), and in that process, one has used the momentum-space version of the BK equation, Eq. (2.13). In other terms, the "subtracted" version, Eq. (2.14), involves a fine cancellation between two large contributions, and this cancellation works only so long as the dipole S matrix obeys the LO BK equation. Any approximation or numerical error in the solution to the latter may spoil this cancellation and thus wash out the equivalence with the original, "unsubtracted", expression (2.11). Indeed, the numerical calculations in [1] have confirmed the equivalence of the two forms in a wide range of transverse momenta. But at the same time, they have shown that, due to numerical errors, some small oscillations persist for high transverse momenta when using the subtracted form. This means that for practical purposes, one is guaranteed to get more stable results when using the unsubtracted form of the quark multiplicity.
Neither Eq. (2.11) nor Eq. (2.14) correspond to the standard k T factorization, since the dipole S matrices are evaluated at the "floating" scale XðξÞ given in Eq. (2.10). In order to arrive at the k T factorized formula presented in [21,22], which is local in X, certain approximations need to be made. First, one observes that due to the subtraction in Eq. (2.14), the integral is dominated by small values ξ ≪ 1. Hence, to the NLO accuracy of interest, it is justified to (i) replace the rapidity argument of the dipole S matrix by its value at ξ ¼ 0, i.e., Sðk; XðξÞÞ → Sðk; X g Þ, and (ii) neglect X g =X 0 ≪ 1 in the upper limit of the integral over ξ. These approximations yield where all the S matrices (explicit or implicit) in the rhs are now evaluated at the rapidity scale X g , that is, for the same rapidity separation from the target as that of the original quark. Equation (2.15) is not any more equivalent to Eqs. (2.11) and (2.14) and, despite the seemingly reasonable approximations, it is rather pathological as it rapidly becomes negative when increasing the transverse momentum of the produced quark (see, e.g., the numerical results in [1]). The reason is that the replacement XðξÞ → X g in the argument of the dipole S matrix leads to an oversubtraction: the negative contribution proportional to Kðk; ξ ¼ 1; X g Þ becomes too large in magnitude and overcompensates for the LO piece in Sðk; X g Þ. Moreover, the extension of the upper limit from 1 − X g =X 0 to 1 is not physically motivated, since it violates constraints imposed from the correct kinematics, and thus, it contains spurious contributions.

C. The running coupling prescription:
Why is this a problem As explained in the Introduction, the experience with the BFKL and BK equations demonstrates that the running coupling (RC) corrections are large for the high-energy evolution, to the point that they should better be viewed as a part of the LO formalism. There are several reasons for that.
The nonlocality of the BFKL/dipole kernels, which is further enhanced by the BFKL diffusion, implies that widely separated transverse (coordinate or momentum) scales, with a priori different values for the running coupling, can contribute to the evolution at a same, given, scale. Furthermore, the effects of the running are amplified by the evolution. For instance, the saturation exponent λ s ≡ d ln Q 2 s =dY, which is a main prediction of the BK equation, is roughly 2 to 3 times smaller when computed with a running coupling than with a fixed coupling. Accordingly, the evolutions of Q s , hence of the gluon distribution, in the two scenarios-fixed coupling vs running coupling-deviate exponentially from each other with increasing Y. For such reasons, it is crucial to use the RC version of the BK equation [11,[37][38][39] when computing the quark multiplicity, both at leading order, cf. Eq. (2.2), and at next-to-leading order, cf. Eq. (2.7). This however introduces complications that we shall discuss in this and the next coming sections.
For consistency with the evolution equation, the explicit factorᾱ s appearing in Eqs. (2.8)-(2.9) for J and J v , which controls the emission of the primary gluon, must be treated as a running coupling too. Since Eqs. (2.8)-(2.9) are explicitly written in transverse momentum space, it seems natural that the scale for the running of that coupling must be a suitable combination of the transverse momenta involved in the emission vertex. When the momentum k ⊥ of the outgoing quark is sufficiently hard, k ⊥ ≳ Q s ðX g Þ, one can simply choseᾱ s ðk ⊥ Þ. This choice is also consistent with other kinematical approximations underlying Eqs. (2.8)-(2.9), notably the relation (2.10) for XðξÞ (see [30] for details). Conversely, the BK equation is generally formulated and solved in transverse coordinate space, cf. Eq. (2.6), where it is more natural to choose a scale for the RC which is a combination of the dipole sizes involved in the splitting vertex (see below for some examples). Albeit perhaps unaesthetic, such a mixed choice for the RC prescriptions-a momentum-space prescription in the calculation of J and J v , but a coordinate-space prescription for the high-energy evolution-is not necessarily a problem: as we shall see, it still allows for meaningful results. However, this has the drawback to introduce some conceptual ambiguities, which ultimately affect the predictive power of the NLO calculation. Specifically, the following issues have been identified by previous studies [1,30]. (b) By the same argument as above, the equivalence between the "unsubtracted" and "subtracted" versions of the NLO result for the multiplicity, cf. Eqs. (2.11) and (2.14), is violated when choosing different RC prescriptions for the primary vertex and the BK equation, respectively. In view of the fine-tuning inherent to the construction of the "subtracted" formula (2.14), this mismatch may have important consequences. As argued in [30], it may contribute to the "negativity problem" which afflicts the original calculation by CXY. Once again, this is confirmed by the numerical study in Ref. [1]. The calculation using the "unsubtracted" formula together with a mixed RC prescription leads to results which are stable and physically meaningful: as in the fixed coupling case, the quark multiplicity receives a negative correction at NLO level, but it remains positive and smooth for all final momenta k ⊥ . On the contrary, the "subtracted" formula leads to a pathological result, which suddenly turns negative at some intermediate value of k ⊥ . This will be further discussed in Sec. III [see, in particular, Fig. 6(b)].
These considerations motivated the recent proposal in Ref. [1] (see the Appendix there) to reformulate the calculation of the NLO multiplicity fully in coordinate space. This should make it possible to use RC prescriptions, which are consistent between the BK equation and the NLO impact factor, and thus remove the ambiguities mentioned at points (a) and (b) above. Moreover, this would allow one to make contact with the phenomenology of DIS at HERA (e.g., in order to constrain the initial condition for the BK equation). Such a reformulation is indeed possible, at least for the N c terms under consideration (as we shall see, the situation is somewhat different for the C F terms): indeed, Eqs. (2.8)-(2.9) can be identically rewritten as the following Fourier transforms:  where the couplingᾱ s has been inserted inside the double integral over r and x, to emphasize that, when using a coordinate-space prescription, this coupling may depend upon all the available transverse-coordinate scales and upon ξ. Some representative diagrams are shown in Fig. 1  Sðr; X g Þ ¼ Sðr; X 0 Þ þ 2πᾱ s ðr; xÞ r 2 x 2 ðr þ xÞ 2 ½Sð−x; XÞSðr þ x; XÞ − Sðr; XÞ where the rewriting in the second line emphasizes that this is the same as the Fourier transform of Eq. (2.13).
The explicit calculations of the one-loop RC corrections to the BK equation [11,[37][38][39] show that these corrections are potentially large whenever there is a strong disparity between the sizes of the three dipoles (the parent and the two daughter ones) involved in the splitting. One can minimize these corrections by choosing the argument of the running coupling as the size of the smallest dipole: α s ðr; xÞ ¼ᾱ s ðr min Þ, with r min ≡ minfjrj; jxj; jr − xjg. But when the dipole sizes are comparable with each other, there is some arbitrariness which is formally an NNLO effect. Besides the aforementioned smallest dipole prescription, one can also use other prescriptions, like the Balitsky prescription [39] or the fastest apparent convergence [19]. Still, one should be rather careful with some choices: for example, for certain values of the dipole sizes the Balitsky prescription takes unphysically small values, and one may question if the contribution from the respective phase space is properly computed (see the discussion in [19]).
The coordinate-space prescriptions used in Eqs. (2.16)-(2.17) for the NLO impact factor and, respectively, the BK equation (2.18) become consistent with each other if the respective running couplings coincide with each other in the limit ξ → 1. When this condition is satisfied, one recovers, e.g., the equivalence between the "unsubtracted" and "subtracted" versions of the quark multiplicity at NLO. This was indeed the case for the RC prescriptions (the Balitsky prescription and a certain generalization of it for ξ ≠ 1, whose precise form is not important for our purposes) used in the coordinate-space calculation presented in the Appendix of [1]. However, the results of that particular calculation turned out to be extremely surprising [and in particular very different from those obtained within the same paper by using the momentum-space RCᾱ s ðk ⊥ Þ within the "unsubtracted" scheme]: the NLO corrections not only change sign (i.e., they become positive), but they are also unacceptably large (by 1 to 2 orders of magnitude larger than the LO result, depending upon the value of k ⊥ ). Since this problem looks identical for both the "unsubtracted" and the "subtracted" formulations of the NLO result, it is clear that its origin must be different from that of the "negativity" problem previously discussed. In the next sections, we shall clarify the origin of this problem and also propose a new RC prescription in coordinate space which avoids that problem and leads to physically meaningful results. Yet, as we shall also explain, that particular prescription is still not entirely satisfactory and cannot be viewed as the ultimate solution to the problem of writing the cross section fully in coordinate space.

III. COORDINATE-SPACE PRESCRIPTIONS FOR THE RUNNING COUPLING
In this section, we shall explore in more detail the coordinate-space calculation of the NLO quark multiplicity based on Eqs. (2.16)-(2.18), with the purpose of clarifying the puzzle raised by the results in the Appendix of [1]. The problem becomes more severe-in the sense that the deviation from the expected physical results (as obtained with a momentum-space running coupling) becomes larger and larger-when increasing k ⊥ (see Fig. 6 in Ref. [1] and also Fig. 3 below). For that reason, in what follows, we shall concentrate on the high-momentum tail of the quark distribution, at k ⊥ ≫ Q s (but our arguments will marginally hold down to k ⊥ ∼ Q s ). For such hard momenta, the exact form of the running coupling used in [1] is not needed, and it just suffices to realize that it reduces tō α s ðr ⊥ Þ. At a first glance, r ⊥ looks as a very reasonable length scale for the argument of the coupling: indeed, for k ⊥ ≫ Q s , the momentum space choice is unambiguouslȳ α s ðk ⊥ Þ, and r is the variable conjugate to k. However, as we show in the following, the running coupling choice and the Fourier transform do not "commute", and therefore the choiceᾱ s ðr ⊥ Þ turns out to be inappropriate.

A. A fake potential
To illustrate the issue, let us start by exploring a simpler example, which is still very close to our actual physical problem 8 Specifically, we consider the following two quantities: Note that both quantities are functions only of k, and the index just refers to the form of the coupling being used. We will evaluate the above expressions using the MV model as an input, i.e., with 8 In fact, it is very similar to one of the C F terms discussed in Sec. IV. ð3:3Þ where Q s is essentially 9 the target saturation momentum and Λ is the usual QCD scale. It is implicitly assumed that the coordinate-space coupling in Eq. (3.2) has been regularized in such a way to avoid the Landau pole at r ⊥ ¼ 1=Λ and that the corresponding regularization is smooth enough to avoid artifacts (like oscillations) when performing the Fourier transform. In practice, whenever an explicit regularization is needed, we will usē whereb ¼ ð11N c − 2N f Þ=12N c and C a positive constant (a free parameter to be varied in the numerical simulations) that was introduced to ensure a smooth and regular behavior of the RC in the "infrared" (i.e., for large r ⊥ ). We consider the regime k ⊥ ≫ Q s , and we first focus on the quantity N k in Eq. (3.1). The result will be dominated by the single scattering term, since multiple scattering contributions will be suppressed by extra powers of Q 2 s =k 2 ⊥ . Notice that the zeroth order term in the expansion of the exponential in Eq. (3.3), i.e., the unity (corresponding to no scattering), leads to a term proportional to δðk ⊥ Þ, which of course vanishes for our purposes. Keeping only the first order term in the expansion (i.e., the contribution of a single scattering), we can write This is the expected result for the high-k ⊥ tail of the distribution of a produced quark within the MV model, including the factorᾱ s ðk ⊥ Þ introduced by hand in the definition (3.1). Now we move on to the hard-k ⊥ regime for the quantity N r in Eq. (3.2), and, similarly to the above, one can check that the single scattering contribution is again proportional to 1=k 4 ⊥ . However, this is not the dominant term anymore, since, due to the inhomogeneity in r ⊥ of the coupling, the unit piece of the S matrix yields a nonzero result. Keeping only that unit piece, our problem amounts to evaluating the Fourier transform of the running coupling, that is ð3:7Þ To write the second line above, we have used the fact that when k ⊥ is very large, so is 1=r ⊥ ∼ k ⊥ , and we can neglect the constant C next to the logarithm. Strictly speaking, the above integration should be restricted to Λ=k ⊥ ≪ k ⊥ r ⊥ ≪ k ⊥ =Λ for the expansion to be valid. In practice though, the lower limit can be safely set to 0 since Λ=k ⊥ ≪ 1, and the integrand does not have a strong support when r ⊥ → 0, so the error will be power suppressed. Furthermore, the upper limit can be set to infinity because now the integration is converging despite the fact that we have ignored the regularization constant C. We thus find ð3:8Þ and we immediately notice that N r , when compared to the correct expression N k given in Eq. which emerges in the absence of an actual scattering and moreover dominates over the respective tail ∼1=k 4 ⊥ introduced by a single scattering (here treated in the MV model). In evaluating N k , the physical picture is indeed correct, as the transverse momentum k ⊥ has been provided (via a 2gluon exchange) by the target potential. On the contrary, when evaluating N r as shown in Eqs. (3.6)-(3.8), a hard momentum k ⊥ emerges due to the singular behavior of α s ðr ⊥ Þ when r ⊥ → 0. In other words, the running coupling is acting like a "fictitious potential," and this cannot be the correct physical picture.

B. From parent dipole to daughter dipole prescriptions
We now move to the actual problem under consideration, i.e., the Fourier transforms which appear in the NLO calculation of the quark multiplicity. As in the previous sections, we restrict to the N c terms for the time being and postpone the discussion of the C F terms to the next section. We therefore need to evaluate J ðk; ξÞ and J v ðk; ξÞ given by Eqs. (2.16) and (2.17), respectively, and we will focus on the real term. For convenience, we will omit the rapidity arguments of J , J v , and S.
Before we discuss various RC prescriptions, let us establish what is the range in the transverse coordinate of the primary gluon, x, which controls the integrations in (2.16). Recall that we consider the production of a relatively hard quark, with a transverse momentum k ⊥ ≫ Q s . Then the Fourier transform selects a relatively small value r ⊥ ∼ 1=k ⊥ ≪ 1=Q s for the parent dipole size. It is then easy to check that the dominant contribution to the integral over x at a fixed (and small) value of r comes from relatively large daughter dipoles, with transverse sizes x ⊥ ≃ jx þ rj ≫ r ⊥ . This is due to the behavior of the kernel in the respective integral, which vanishes as x ⊥ → 0 and decays only very slowly, as 1=x 2 ⊥ , when x ⊥ → ∞. Of course, the integral will be eventually cutoff by the physics of saturation, i.e., by the fact that any of the S matrices appearing in Eq. (2.16) will rapidly vanish when increasing x ⊥ above 1=Q s . But the dominant contribution comes from the range r ⊥ ≪ x ⊥ ≲ 1=Q s and more precisely. from its upper limit, i.e., from x ⊥ ∼ 1=Q s , due to the fact that the scattering amplitude Tðx ⊥ Þ ≡ 1 − Sðx ⊥ Þ grows quadratically with x ⊥ so long as x ⊥ ≪ 1=Q s . We thus conclude that the would-be dominant contribution to the integral over x, as coming from very large daughter dipoles, is independent of r, and hence it does not at all contribute to the subsequent Fourier transform, The latter is rather dominated by daughter dipoles whose sizes are comparable to that of their parent: x ⊥ ∼ r ⊥ . For such small dipoles, one can use the single scattering approximation within the MV model and thus deduce which is indeed a correction of Oðᾱ s Þ to the tree-level result at high k ⊥ , cf. Eq. (3.5). A similar estimate holds for the virtual term J v ðk; ξÞ. Notice that some of the previous arguments have exploited the fact that ξ ≠ 1, as is indeed the case for the generic values contributing to the integral over x ¼ 1 − ξ in Eq. (2.7). The limit ξ → 1 of these results is well defined and is discussed below. To summarize, the double integral in Eq. (2.16) is not controlled by the typical gluon emissions, which are favored by the emission kernel and are relatively soft (p ⊥ ∼ Q s ), but by the rare emission of a hard gluon with p ⊥ ∼ k ⊥ . This conclusion could have been anticipated on physical grounds: the relatively large momentum k ⊥ ≫ Q s of the measured quark cannot be provided by its scattering off the nuclear target; hence, it must be balanced by the emission of an equally hard primary gluon, with transverse momentum p ⊥ ≃ k ⊥ . In coordinate space, this corresponds indeed to a transverse size x ⊥ ∼ r ⊥ . This correct physical picture should be preserved by any prescription used for the running of the coupling.
First of all, it is clear that the whole previous discussion remains unchanged if in Eqs. (2.16) and (2.17), one uses a running coupling with the scale set by the "external" momentum k ⊥ . This choiceᾱ s ðk ⊥ Þ is physically well motivated when k ⊥ ≫ Q s , as already explained. So, in what follows, we shall use the results computed with the momentum-space RCᾱ s ðk ⊥ Þ as a benchmark for analyzing the suitability of other RC prescriptions formulated in coordinate space. For clarity, we denote the former with a superscript "mom" and the latter with a superscript "pos".
Since the Fourier transform in Eq. (2.16) selects parent dipole sizes r ⊥ ∼ 1=k ⊥ , it might seem reasonable to use the RCᾱ s ðr ⊥ Þ. This was essentially the choice made in [1] (at least for sufficiently high k ⊥ ). Yet, as shown by the numerical results presented in the Appendix of Ref. [1], this prescription leads to physically unacceptable results. The problem is further illustrated in Figs. 3(a) and 3(b). Specifically, Fig. 3(a) shows the ratio J pos =J mom between the results for the "real" piece J obtained with the coordinate-space prescription 10ᾱ s ðr ⊥ Þ and the momentum-space prescriptionᾱ s ðk ⊥ Þ ≡ 1=½b lnðk 2 ⊥ =Λ 2 Þ. For the dipole S matrix, we have used the solution to rcBK with smallest-dipole prescriptionᾱ s ðr min Þ and with an initial condition given by the MV model, evolved from X 0 ¼ 10 −2 down to X ¼ 10 −3 . When evaluating J ðk; ξ; XÞ according Eq. (2.16), we have treated X, ξ, and k ⊥ as independent variables. 11 Similarly, Fig. 3(b) shows the corresponding ratio J pos v =J mom v for the "virtual" piece J v . By inspection of these figures, it is manifest that the results corresponding to the two RC prescriptions are inconsistent. For instance, instead of being close to one for sufficiently high k ⊥ , the ratio J pos =J mom is negative, and its modulus is very large (≫ 1) and rapidly increasing with k ⊥ , according to a power law.
These results demonstrate the failure of the coordinatespace prescriptionᾱ s ðr ⊥ Þ for the NLO impact factor. This failure is easy to understand in light of our discussion in this and the previous sections. Recall that, by itself, the kernel in Eq. (2.16) would favor large daughter dipoles with x ⊥ ≫ r ⊥ (meaning soft primary emissions), which however are eliminated by the final Fourier transform r → k, since their contribution is quasi-independent of r, cf. (3.9). However, the situation changes in the presence of a RC likeᾱ s ðr ⊥ Þ, which explicitly depends upon r ⊥ : this coupling effectively acts as a "potential," which is singular as r ⊥ → 0, by asymptotic freedom, meaning that it can transfer an arbitrarily hard transverse momentum to the incoming quark. With theᾱ s ðr ⊥ Þ prescription, the Fourier transform hence introduces a spurious tail ∝ 1=k 2 ⊥ with a wrong sign, cf. Eq. (3.8), which at large k ⊥ dominates over the physical tail ∝ 1=k 4 ⊥ -in agreement with the numerical results exhibited in Figs.

3(a) and 3(b).
At this stage, one may wonder about the implications of our previous discussion for the solution to rcBK, where it is quite common to use a RC prescription-likeᾱ s ðr min Þ or the Balitsky prescription discussed in Appendix A-which becomes equivalent to the parent-dipole prescription α s ðr ⊥ Þ whenever r ⊥ is the smallest scale. In that context, such a prescription is well-known to be well-behaved. We show now that this is not in contradiction with our above arguments that J and J v are ill-behaved when computed with theᾱ s ðr ⊥ Þ prescription. The crucial point is that the rhs of the BK equation (2.18) involves the difference J ðk; ξ; XÞ − J v ðk; ξ; XÞ with ξ ¼ 1, that is where we have directly inserted a RCᾱ s ðr ⊥ Þ and we have omitted the rapidity arguments for simplicity. The main difference compared to the previous discussion is that now the unphysical contributions due to large dipoles cancel out in this difference. Indeed, Eq. (3.11) involves the dipole kernel, which decays much faster-namely, as r 2 ⊥ =x 4 ⊥ -at large x ⊥ ≫ r ⊥ . Using again the MV model for the dipole S matrix and restricting ourselves to the single scattering approximation, one finds  The actual dependence of X on ξ and k ⊥ , shown in Eq. (2.10), is unessential for the present argument.
where we have kept only the contribution from relatively large daughter dipoles, within the range r ⊥ ≪ x ⊥ ≪ 1=Q s . (This is indeed the dominant contribution, since it is enhanced by an additional transverse logarithm.) The integrand in Eq. (3.12) is very different from that occurring in Eq. (3.9) for the "real" piece J alone. Compared to Eq. (3.9), the integration over x now introduces an extra logarithmic dependence on the parent dipole size r ⊥ , on top of that encoded in the RCᾱ s ðr ⊥ Þ. This additional dependence is the actual physical source for the transverse momentum k of the produced quark. The final Fourier transform in Eq. (3.12) is in fact similar to that appearing in Eq. (3.10) and thus yields a high-momentum tail ∼1=k 4 ⊥ , as expected for the first step in the BK evolution of the quark multiplicity.
Albeit the previous argument was directly constructed for the limit ξ → 1, as relevant for the BK equation, it actually holds for generic values of ξ: the difference J ðk; ξ; XÞ − J v ðk; ξ; XÞ is free of the "fake potential" problem for any value ξ ≤ 1. This can be easily verified by making a change of variables in Eq. (2.17) in order to match the arguments of the various S matrices there with those in Eq. (3.11). This is further confirmed by the numerical results in Fig. 3(c). This observation also shows that the ultimate reason why the calculation of the NLO multiplicity using the parent dipole prescriptionᾱ s ðr ⊥ Þ is so ill behaved is because the "real" and "virtual" contributions to Eq. (2.2) are differently weighted by the quark distribution, which is taken at a longitudinal momentum fraction equal to x p =ξ in the "real" contribution, but equal to x p in the "virtual" one.
The previous discussion immediately suggests a better RC prescription to be used in coordinate space: the scale of the RC should rather be set by the daughter dipole size x ⊥ . Indeed, with the choiceᾱ s ðx ⊥ Þ, the spurious contribution generated by large (x ⊥ ≫ r ⊥ ) daughter dipoles remains independent of r, meaning it is eventually removed by the Fourier transform, as it should. Then the net contribution comes solely from the physical configurations with x ⊥ ∼ r ⊥ . For such configurations, the coordinate-space prescription α s ðx ⊥ Þ should be equivalent to the momentum-space oneᾱ s ðk ⊥ Þ. This is confirmed by the numerical results shown in Figs. 4(a) and 4(b), which are the counterpart of those in Figs. 3(a) and 3(b), but with a daughter-dipole prescription α s ðx ⊥ Þ in coordinate space. As visible in these figures, the results are now nicely consistent between the coordinatespace and the momentum-space prescriptions. The remaining difference should be viewed as a measure of the RC scheme dependence of our calculation, and this dependence turns out to be quite small (and decreasing with increasing k ⊥ , as expected). Similarly, the behavior of J − J v , shown in Fig. 4(c) for the daughter-dipole prescriptionᾱ s ðx ⊥ Þ, is also well-behaved and within our scheme dependence.

C. NLO results with various running coupling prescriptions
To further illustrate our discussion in this section, we exhibit in Figs. 5 and 6 numerical results for the NLO multiplicity obtained with different prescriptions for the RC. For more clarity, the NLO results are normalized by the respective LO predictions, as obtained by using Eq. (2.2) together with the solution to the BK equation with running coupling. The RC prescription for rcBK is the minimal dipole prescriptionᾱ s ðr min Þ, unless otherwise specified. As before, we use the MV model as an initial condition at X 0 ¼ 10 −2 and present the results for the cross section as a function of k ⊥ corresponding to a COM energy ffiffi ffi s p ¼ 500 GeV and a pseudorapidity η ¼ 3.2 for the produced quark. The corresponding values for x p and X g can be read off Eq. (2.5). All the curves presented in Fig. 5 are obtained by using the "unsubtracted" formula Eq. (2.11). The four curves in Fig. 5(a) correspond to four different choices for the coupling, which is explicit in the NLO impact factor (i.e., the coupling which controls the emission of the primary gluon and whose running has been the main issue of this section): (i) a fixed couplingᾱ s ¼ 0.2, (ii) a one-loop running coupling with a momentum-space prescription, α s ðk ⊥ Þ, (iii) a one-loop RC with a parent dipole prescription,ᾱ s ðr ⊥ Þ, and finally (iv) a one-loop RC with a daughter dipole prescription,ᾱ s ðx ⊥ Þ. The regulariation used for the one-loop RC is as shown in Eq.  Fig. 3 except that the coordinate-space prescription used for the RC coupling which controls the primary emission is the daughter dipole prescriptionᾱ s ðx ⊥ Þ. All the ratios are now rather close to one (at least for sufficiently large values of k ⊥ , where the present approximations are supposed to apply).
chosen in such a way thatᾱ s ðr ⊥ Þ → 0.67 as r ⊥ → ∞. We observe that the two curves corresponding toᾱ s ðk ⊥ Þ and tō α s ðx ⊥ Þ are not only very close to each other (thus confirming the little scheme dependence in our calculation), but also very close to the corresponding results at a fixed coupling. That is, the effect of the running of the coupling is not so important in so far as the primary emission is concerned. This is due to the fact that the running of the coupling is merely logarithmic and the interval in k ⊥ that we consider in Fig. 5 is quite small. Besides, the effects of this particular coupling are not amplified by the evolution, unlike what happens for the coupling occurring in the BK equation. On the other hand, the results obtained with the parent-dipole prescription α s ðr ⊥ Þ are dramatically different, in agreement with the previous discussion in this section (and also with the original results in [1]): they differ from the correct results by up to 2 orders of magnitude, and this difference keeps increasing with k ⊥ , due to the spurious high-k ⊥ tail introduced by the Fourier transform ofᾱ s ðr ⊥ Þ. The curve labeled as "ᾱ s ðk ⊥ Þ" in Fig. 5(b) is exactly the same as the respective curve in Fig. 5(a). But the other curve in Fig. 5(b), denoted as "ᾱ s ðx ⊥ Þ," is now obtained by systematically using the daughter dipole prescription α s ðx ⊥ Þ throughout the calculation, including within rcBK. That is, it differs from the respective curve in Fig. 5(a) by the RC prescription used for the dipole evolution. The two curves in Fig. 5(b) are quite close to each other, albeit their difference is somewhat larger than that observed between the three "physical" curves in Fig. 5(a) [notice however the different vertical scales used in these two figures]. Such a difference was in fact to be expected: a RC with a daughter dipole prescription gives a stronger weighting to the emission of gluons with a soft transverse momenta (i.e., to large daughter dipoles) and thus leads to a somewhat faster evolution. We therefore expect the difference between the two curves in Fig. 5(b) to grow with increasing the phase space for energy evolution.
Having a unique RC prescription in the whole NLO calculation looks conceptually appealing, in that it renders the calculation more homogeneous. In particular, this removes both ambiguities mentioned at points (a) and (b) in Sec. II C, as numerically demonstrated in Fig. 6. The curves "ᾱ s ðk ⊥ Þ" and "ᾱ s ðx ⊥ Þ" shown in this figure are obtained by using the same RC prescriptions as in Fig. 5(b). As expected, the curve "ᾱ s ðx ⊥ Þ," which uses a unique, daughter-dipole prescription in all the stages of the calculation, has an unambiguous LO limit [cf. Fig. 6(a)] and yields identical results with both the "unsubtracted" and the "subtracted" formulations of the NLO multiplicity [cf. Fig. 6(b)]-at variance with the curve labeled "ᾱ s ðk ⊥ Þ", which uses a mixed set of RC prescriptions.
This being said, the calculation underlying the curve "ᾱ s ðk ⊥ Þ" in Fig. 5(b), which uses a mixed RC prescription together with the "unsubtracted" formula (2.11) for the cross section, is still the one to be a priori trusted at high energy. This is so since, strictly speaking, the daughter dipole prescriptionᾱ s ðx ⊥ Þ is not fully suitable for the BK equation, in that it artificially accelerates the evolution. To understand that, recall that the dominant dipole configurations are not the same for the calculation of the NLO impact factor (i.e., for the emission of the primary gluon) and for the BK evolution (i.e., for the subsequent emissions of soft gluons, as resummed by the BK equation). Indeed, unlike the primary gluon, which is as hard as the produced quark (so the associated "daughter dipole" has a small size x ⊥ ∼ r ⊥ , as previously discussed), the subsequent emissions in the "hard-to-soft" BK evolution involve predominantly gluons with smaller and smaller transverse momenta, or "daughter dipoles" with larger and larger transverse sizes. For such typical emissions, pQCD (and in particular the experience with the DGLAP equation [40]) instructs us that the proper scale to be used in the RC is the transverse size of the parent dipole (see also the discussion in Appendix A below). By asymptotic freedom, this typically leads to a smaller value for the coupling as compared to the daughter dipole prescription, and hence to a somewhat slower evolution. A further argument in favor of the momentum-space prescriptionᾱ s ðk ⊥ Þ will be provided by the discussion of the C F terms in the next section.

IV. THE C F TERMS
In this section, we shall extend the previous discussion to the remaining contributions to the quark multiplicity at NLO, those proportional to the quark Casimir C F . If we preferred to discuss these terms separately, it is not because they were exempted of the "fake potential" problem aforementioned-as we shall see, they suffer from a similar problem when computed in coordinate space and with a parent dipole running-coupling prescription-but rather because they encode a different physical regime. Unlike the N c terms, which include the small-x limit of the primary gluon emission and thus overlap with the LO BK evolution, these C F terms include the respective collinear limits-i.e., the limits where the primary gluon is collinear with either the incoming quark, or with the outgoing one-and hence overlap with the DGLAP evolution of the quark distribution function, and of the quark-into-hadron fragmentation function, respectively. These collinear limits signal themselves via (infrared) logarithmic divergences, which need to be regularized and subtracted away (since the respective contributions have already been included via the above mentioned DGLAP evolution). In practice, this is conveniently done by using dimensional regularization together with minimal subtraction. The remaining, finite, terms are pure NLO corrections to the impact factor. This renormalization procedure is explained in detail in Ref. [21,22], from which we shall simply quote the necessary results.
The contribution of the C F terms to the quark multiplicity at NLO can be written as ð4:1Þ where I and I v are "real" and "virtual" terms associated with the same process as described before: the noneikonal splitting q → qg followed by the scattering off the nucleus. For more clarity, we first present the original version of these terms, prior to renormalization. When computed in transverse momentum space, they read The collinear divergences (at q ¼ k and q ¼ k=ξ in the "real" term and, respectively, at q ¼ k and q ¼ ξk in the "virtual" term) are indeed manifest. It is a straightforward exercise in Fourier transforms to write the above formulas as double integrations in coordinate space, so that they look more similar in structure to those in Eqs. (2.16) and (2.17). For example, for the "real" term we find  x · ðx þ rÞ x 2 ðx þ rÞ 2 ½Sðr; XðξÞÞ þ Sðξr; XðξÞÞ − 2Sðξr − ð1 − ξÞx; XðξÞÞ; ð4:4Þ where the three terms within the square brackets correspond to the terms obtained by expanding the square in the integrand of Eq. (4.2) and also to the three diagrams shown in Fig. 7. The first two diagrams, Figs. 7(a) and 7(b), are the expected counterparts of the two diagrams shown in Fig. 1, corresponding to the other possible time orderings for the emission of the primary gluon. In the eikonal limit ξ → 1, these four diagrams would together generate the first step in the BK evolution. This last statement also shows that, by themselves, the first two terms within the brackets in Eq. (4.4)-corresponding to the squared terms, 1=ðk − qÞ 2 and 1=ðk − ξqÞ 2 , in Eq. (4.2)-would also generate a longitudinal logarithm via the ξ → 1 limit of the integral in Eq. (4.1); that is, they would contribute to the high-energy evolution as well. The third term in Eq. (4.4)-corresponding to the "crossed" term in the square in Eq. (4.2)-has been introduced by hand in Refs. [21,22] in order to subtract the ξ → 1 limit of the first two terms. More precisely, this term has been added to the N c terms and subtracted from the C F terms, in such a way to disentangle the soft from the collinear divergences. The apparent mismatch between the color factors of the "added" and "subtracted" terms, N c and, respectively, C F =2, is irrelevant in the large N c limit of interest here. But even for finite N c , this mismatch is precisely compensated by that piece of the diagram in Fig. 1(a) which would be suppressed at large N c . This piece is illustrated in Fig. 7(c) (see the caption to Fig. 7 for more details). A similar discussion applies to the "virtual" term I v .
To summarize, the "real" and "virtual" C F terms, I and I v , separately vanish in the eikonal limit ξ → 1, thus guaranteeing that the integration in Eq. (4.1) will not generate any large longitudinal logarithm. They also vanish-once again, due to the subtraction of the third term in Eq. (4.4)-in the absence of any scattering, i.e., when SðrÞ ¼ 1 or, equivalently, SðkÞ ¼ ð2πÞ 2 δ ð2Þ ðkÞ. A similar property holds for the N c terms, as manifest, e.g., in Eqs. (2.16)-(2.17). This property is important since it guarantees that there is no particle production in the absence of scattering, as expected on physical grounds.
However, Eqs. (4.2), (4.3), and (4.4) are afflicted with infrared divergences associated with collinear emissions. In momentum space, these divergences were already noticed after Eqs. (4.2)-(4.3). Importantly, they refer to the squared terms alone (in either I and I v ); the respective crossed terms are free of singularities for any ξ < 1. In coordinate space, e.g., in Eq. (4.4), the infinities occur in the limit where the gluon is emitted very far from the quark, that is when x ⊥ gets very large. Once again, this happens only for the first two terms in the square bracket in Eq. (4.4), corresponding to the first two diagrams in Fig. 7. In Refs. [21,22], the subtraction of the collinear divergences has been performed by working in momentum space, i.e., at the level of Eqs. (4.2)-(4.3), and for a fixed coupling. The final results after subtraction can be written as . 7. Feynman graphs contributing to the C F term I in coordinate space, which illustrate the three terms in the brackets in Eq. (4.4). In graph (a), the interactions between the gluon and the shockwave fully cancel between the direct amplitude and the complex conjugate amplitude, by unitarity: Vð−xÞV † ð−xÞ ¼ 1 (here, the Wilson line is in the adjoint representation). In graph (b), the gluon does not interact with the shockwave neither in the amplitude nor in the complex conjugate one. The "photonlike" line replacing the primary gluon in graph (c) represents the color singlet piece of that gluon (hence it does not interact with the shockwave). Indeed, graph (c) represents only a piece of an actual Feynman graph, namely that piece of the diagram in Fig. 1(a) which is suppressed at large N c (the singlet piece in the color decomposition ½N 2 c − 1 ¼ ½N c × N c − ½1 of the gluon). Following Refs. [21,22], the overall factor 1=2N c multiplying this suppressed piece has been conveniently rewritten as 1=2N c ¼ N c =2 − C F , i.e., as the difference between two unsuppressed contributions. This rewritting generates two terms, one proportional to N c =2 that was included in the N c terms [as the first term in the rhs of Eq. (2.8) for J ], the other one proportional to C F , which is the third term in Eq. (4.4  with μ a factorization scale and Besides being finite, the above expressions for I fin and I fin v preserve the "good" properties of their unregularized versions I and I v that we previously discussed. First, they vanish in the absence of any scattering and for k ≠ 0; this is obvious for I fin v and can also be checked for I fin . [When SðrÞ ¼ 1 and SðqÞ ¼ ð2πÞ 2 δ ð2Þ ðqÞ, we can easily perform the integrations in Eq. (4.5) and thus find that I fin ¼ 0.] Second, in Appendix B, we check that when using the finite expressions for the C F terms in Eqs. (4.5) and (4.6) to compute the integral over ξ in Eq. (4.1), one does not generate any longitudinal logarithm. As shown by the above arguments and also by the manipulations in Appendix B, these "good" properties are guaranteed to hold with either a fixed coupling, or a running couplinḡ α s ðk ⊥ Þ with the scale set by the external transverse momentum. However, as we shall shortly argue, they would not hold anymore if the scale for the RC is rather chosen in coordinate space.
Let us first observe that a coordinate-space RC prescription would not be very natural in the present context, where the subtraction of the collinear divergences has been performed in momentum space. Indeed, the terms affected by this subtraction [the first term in Eq. (4.5) and the whole term in Eq. (4.6)] cannot be written as a double integration in coordinate space-unlike the original, singular, expressions for I and I v , cf. Eq. (4.4). As a consequence, we have no control anymore on the sizes of the "daughter" dipoles [the coordinate x in Eq. (4.4)]. Thus, if one insists in using a running coupling with a transverse coordinate scale, then the only option isᾱ s ðr ⊥ Þ. By inspection of Eqs. (4.5) and (4.6), it is clear that such a choice would entail severe ambiguities. For instance, if one inserts the couplingᾱ s ðr ⊥ Þ inside the integral over r in the second equality of Eq. (4.6), then one violates the equivalence with the expression shown in the first equality there (the a priori result of the collinear subtraction in momentum space). We are facing here exactly the same ambiguity as previously discussed in relation with Eqs. (3.1)-(3.2). From that example, we know that the use of the RCᾱ s ðr ⊥ Þ would introduce the problem of the "fake potential". This problem is illustrated in Fig. 8(b), which shows the ratio I pos v =I mom v between the results for I v obtained with the position-space RCᾱ s ðr ⊥ Þ and the momentum-space RCᾱ s ðk ⊥ Þ; one can see that I pos v has a different sign and also a different power tail at high k ⊥ as compared to I mom v -in agreement with the discussion around Eqs. (3.6)-(3.8).
An equally severe ambiguity, leading again to unphysical results, occurs if one tries to use the position-space RC α s ðr ⊥ Þ in Eq. (4.5) for the "real" term I fin . The last term there is exactly the same as the third term in Eq. (4.4); hence, it can still be written as a double integral over r and x. So, it is in principle possible to use the RCᾱ s ðr ⊥ Þ in all the terms appearing in Eq. (4.5) (inside the respective integrations, of course). However, if one does so, one obtains the results shown in Fig. 8(a), which differ by several orders of magnitude (and also in sign) from those obtained with the momentum-space RCᾱ s ðk ⊥ Þ. Besides, the use of the position-space RCᾱ s ðr ⊥ Þ would spoil the cancellation between the two terms in Eq. (4.5) in the limit where there is no scattering. It would furthermore spoil the  8. (a) The "real" and (b) the "virtual" C F terms, I and I v , computed with the parent dipole prescriptionᾱ s ðr ⊥ Þ and normalized to the corresponding results obtained with the momentum-space prescriptionᾱ s ðk ⊥ Þ. In the middle graph, the values of ξ need not be specified, since the ratio I pos v =I mom v is independent of ξ, as manifest by inspection Eq. (4.6). Also, we exclude the limiting value ξ ¼ 1 from these plots because the densities I and I v are separately divergent in the limit ξ → 1. Moreover, the difference I mom − I mom v in (c) vanishes as ξ → 1, as shown in Appendix B. cancellation between I fin and I fin v in the limit ξ → 1 (thus generating a spurious longitudinal logarithm in the C F terms); see the discussion in Appendix B.
In Fig. 9, we show the complete result for the NLO quark multiplicity as emerging from the present calculations, including both the N c terms and the C F terms, and for three RC prescriptions: fixed couplingᾱ s ¼ 0.2, momentumspace RCᾱ s ðk ⊥ Þ, and the (pathological) coordinate-space RCᾱ s ðr ⊥ Þ. To better appreciate the importance of the C F terms, we also display the NLO result including just the N c terms (with momentum-space RC), that is, the curve that was denoted "ᾱ s ðk ⊥ Þ" in Fig. 5(a). One thus sees that the effect of adding the C F terms is indeed sizeable: the full NLO result is still smaller than the LO one, but the difference between the two is considerably reduced. A similar phenomenon-namely the existence of large C F and N c corrections, which have opposite signs and largely cancel each other-has been also noticed in the recent numerical calculation of DIS at NLO in Ref. [31]. That calculation too was based on a factorization scheme nonlocal in rapidity (similar to the one that we employed in this paper) together with the NLO impact factor for the virtual photon from Refs. [32,33].
To summarize, the problem with a coordinate-space prescription in the RC looks even more severe for the C F terms than for the N c terms discussed previously. Unlike for the latter, we have no control on the size of the daughter dipoles anymore, so one cannot use the couplingᾱ s ðx ⊥ Þ to avoid spurious contributions generated by the Fourier transform of the running coupling. The only reasonable choice seems to be a coupling running with a transverse momentum scale, such as the momentum k ⊥ of the produced particle.

V. CONCLUSION AND OUTLOOK
In this paper, we have continued our work on the NLO corrections to the single inclusive particle production at forward rapidities in high energy proton-nucleus collisions. For simplicity, we considered only the production of a quark, but ignored the gluon channel, as well as the subsequent fragmentation of the quark into hadrons. Our main focus was on the choice of the scale for the running of the QCD coupling. This coupling appears in two distinct places in the calculation: the high-energy evolution of the dipole S matrix, as described by the (NLO) BK equation and the NLO impact factor.
The BK equation is generally solved by using the transverse coordinate space representation, whereas the NLO calculation of the impact factor is most naturally performed in the transverse momentum space. This naturally leads to the use of different running coupling prescriptions-a coordinate-space prescription for the BK equation, but a momentum-space one for the NLO impact factor-which in turn introduces some ambiguities in practice (because the running of the coupling does not "commute" with the Fourier transform). It is then legitimate to ask whether one can eliminate such ambiguities by using a unified RC prescription throughout the calculation. A first proposal in that sense [1], using a coordinate-space formulation for the RC, was found to lead to unphysical results: the NLO correction to the quark multiplicity was larger by up to 2 orders of magnitude than the LO term, and it also had the wrong sign compared to the calculation using the momentum-space RC prescription.
Here, we have been able to understand the origin of the problem with the coordinate-space prescription proposed in [1]. We first observed that in the kinematical regime of interest, namely for relatively large transverse momenta k ⊥ ≳ Q s for the produced quark, that the prescription was essentially the same as the parent-dipole prescription α s ðr ⊥ Þ, with r ⊥ ∼ 1=k ⊥ the transverse size of the parent dipole size. With this prescription, the RCᾱ s ðr ⊥ Þ acts as a spurious potential and leads to a large, but unphysical, contribution to the cross section from primary gluons (daughter dipoles) with large transverse sizes x ⊥ ≫ r ⊥ . This contribution is unphysical in that it violates transverse momentum conservation (i.e., the fact that the momentum of the primary gluon should match the one of the hard outgoing quark). In a correct calculation, this contribution would cancel due to the final Fourier transform, but this cancellation is badly spoiled by the fake "potential"ᾱ s ðr ⊥ Þ, which can transmit a large momentum. Regarding the C F terms, the choiceᾱ s ðr ⊥ Þ leads to additional problems: first, it generates nonzero, and thus unphysical, contributions in the limit of no scattering; second, it spoils the cancellation of the longitudinal logarithms (which in this factorization scheme should exist only in the N c terms).
We have also found that there is a meaningful scale choice in coordinate space, at least so long as the N c terms are considered. This is the "daughter" dipole prescription α s ðx ⊥ Þ, which keeps intact the required cancellations. Numerically, we have seen that this choice gives almost the same results as the choiceᾱ s ðk ⊥ Þ, with the small difference being an acceptable scheme dependence. However, this choiceᾱ s ðx ⊥ Þ is not possible in the C F terms, where the coordinate x ⊥ of the primary gluon has been integrated over prior to the subtraction of the collinear divergencies (and hence this is not present anymore in the finite results). Thus, for these terms, the use of the momentum space couplingᾱ s ðk ⊥ Þ seems to be mandatory. With the above choices, that is withᾱ s ðk ⊥ Þ orᾱ s ðx ⊥ Þ for the N c terms andᾱ s ðk ⊥ Þ for the C F terms, we have numerically computed the quark multiplicity at NLO. We have found that the NLO correction reduces the LO result by 25%-30%, as one expects in a controlled pQCD calculation and in this regime of transverse momenta. Also, we have found that the NLO N c and C F contributions are of the same order, as also observed in the NLO calculation of the structure functions in DIS.