Single inclusive forward hadron production at next-to-leading order

We discuss single inclusive hadron production from a high energy quark scattering off a strong target color field in the Color Glass Condensate formalism. Recent calculations of this process at the next-to-leading order accuracy have led to negative cross sections at large transverse momenta. We identify the origin of this problem as an oversubtraction of the rapidity divergence into the Balitsky-Kovchegov evolution equation for the target. We propose a new way to implement the kinematical restriction on the emitted gluons to overcome this difficulty.


I. INTRODUCTION
Hadronic reactions at modern collider energies reach a kinematical domain where gluon densities can be nonperturbatively large even at short distance scales where the QCD coupling is weak. A convenient effective theory of QCD in this regime is provided by the Color Glass Condensate (see e.g. [1]), which describes the nonlinear small x degrees of freedom in a hadron or nucleus as a classical color field. An ideal way to study these dense color fields is to probe them with simple dilute probes in a high energy collision, such as in deep inelastic scattering or forward particle production in proton-nucleus collisions. In the latter case the dilute probe is provided by the relatively well understood large x partons of the probe proton which, at forward rapidity, scatter off the small x color field of the target.
Several calculations [2][3][4][5][6] of forward single inclusive particle production in this framework have provided a good description of available experimental data using the leading order expression for the cross section [7]. As is often the case in QCD, the leading order calculations leave the overall normalization of the cross section quite uncertain. It would, therefore, be desirable to systematically go to higher orders in perturbation theory in these cross section calculations. Similar developments towards higher order have recently taken place concerning the (Balitsky-Kovchegov (BK) [8,9] or JIMWLK) high energy evolution equations [10][11][12][13][14][15][16][17][18], and DIS cross sections [19,20].
An important advance in pushing the CGC framework to NLO accuracy was the calculation [21,22] of NLO single inclusive particle production in forward protonnucleus collisions (see also the earlier works [2,23,24]). We will here refer to the cross section formulae derived in Refs. [21,22] as the "CXY" result according to the authors. Here it was shown that the divergences in the rapidity (or longitudinal momentum) and transverse momentum integrals appearing in the NLO calculation can be factorized into the BK and DGLAP evolution of the target and projectile, respectively. In a subsequent calculation [25] it was shown that in the large transverse momentum limit the calculation reduces to the appropri-ate tree-level process in collinear factorization, although in this case without a factorization of the rapidity divergence.
In the first numerical implementation [26] of the factorization framework of [21,22] the NLO corrections turned out to be large and negative at large transverse momenta of the produced particles, to the extent that the total cross section becomes negative. A solution to this problematic behavior has been suggested to lie in the detailed implementation of the factorization of the rapidity divergence [27] or in a kinematical constraint that must be imposed on the phase space of emitted gluons [28]. Indeed a recent implementation [29] of this phase space constraint has alleviated the problem, without however removing it completely.
In this paper, we suggest an alternative way of implementing BK-factorization in the calculation of Chirilli et al. [21,22] that combines aspects of these previous works. Our suggestion includes, as in [27], an explicit rapidity factorization scale in the "hard functions" of the NLO part of the cross section. Since the dependence of the cross section on this scale cancels against the rapidity scale to which the target is evolved, the total cross section formulation is explicitly independent of this factorization scale. Here we suggest implementing the "kinematical constraint" or, more precisely, ordering in light cone energy, by making this factorization scale dependent on the transverse momentum of the produced particle. We show that this allows for a renormalization prescription that makes the negative large momentum contribution to the cross section arbitrarily small. Combined with a corresponding form of the BK equation this should in the future make it possible to resum the problematic contributions at large transverse momentum similarly as recently suggested for the NLO BK equation [16,18].
For simplicity we will here only address the quark channel q → q and perform numerical calculations only for the Golec-Biernat and Wüsthoff (GBW) [30] parametrization of the dipole cross section. The gaussian k ⊥ -spectrum at leading order in the GBW model has the advantage for the purpose of this paper of being very clearly distinct from the power-law behavior of the NLO contributions. This exacerbates the problem of negative cross sections arXiv:1604.00225v2 [hep-ph] 3 Jul 2016 and should therefore put any attempt to stabilize the perturbative expansion to a more stringent test. A fuller phenomenological analysis of single inclusive particle production would require implementing also the gluon initiated channel, and a more realistic BK-evolved dipole cross section. Ultimately this should include a solution of the NLO BK equation [16,18] and a simultaneous NLO fit to DIS data. At present we leave these further steps for future work.
This paper is structured as follows. We will first, in Sec. II briefly recall CXY results [21,22] for the NLO corrections to the cross section in the quark-initiated channel and how the factorization into DGLAP and BK evolution is performed there. We will then discuss the introduction of an explicit rapidity factorization scale as advocated in [27] and explicitly show the modification to the NLO spectrum resulting from a variation of this scale in Sec. III. We then, in Sec. IV, discuss imposing the additional "kinematical" constraint of k − -ordering (see e.g. [31]) on the rapidity factorization, which we implement as a momentum dependence of the rapidity factorization scale. Section V then concludes with a brief discussion of possible future steps.

II. SINGLE INCLUSIVE PARTICLE PRODUCTION AT NLO
Our starting point are the CXY formulae derived in Ref. [22]. We will concentrate here on the quark channel, for which we write the CXY result as: Here we have slightly altered the notation of [22] by including transverse momentum integrals in the functions I, J and leaving out an overall integration over the impact parameter b, thus our expression is for the multiplicity and not the cross section. The kinematical variables are defined as p = zk, x p = p ⊥ e y h /(z √ s), τ = zx p , Most important for our discussion here is the momentum fraction ξ: the fragmenting quark carries a fraction ξ of the incoming quark longitudinal momentum. Thus the incoming quark has a momentum fraction x p /ξ of the incoming proton, where x p is the probe momentum fraction in leading order kinematics. The radiated gluon in the NLO terms carries a longitudinal momentum fraction 1 − ξ: i.e. the limit ξ → 1 corresponds to the soft gluon emission that must be resummed into BK evolution of the target. These cross sections are all expressed in terms of the Fourier-transform of the fundamental representation dipole operator The dipole is related to the notation of [22] by an overall integration over the impact parameter In the following we leave out the explicit impact parameter dependence of S(k ⊥ , b) from our notation. Due to the unitarity of the Wilson lines U the dipole cross section satisfies the normalization condition The expressions here use the mean field approximation, replacing the expectation value of the product of two dipole operators by the product of two expectation values S(q ⊥ )S(l ⊥ ). The superscript (0) refers to the fact that at this stage the dipole operator in the leading order part of (1) is the unrenormalized "bare" dipole operator. As noted in [22], several important features are already visible in these expressions. First, the terms with an explicit coefficient C F (i.e. I and I v ) vanish in the limit ξ → 1 but have collinear divergences in the transverse momentum integration. These are nicely treated in [22] by dimensional regularization in the transverse momentum integrals and factorized into the DGLAP evolution of the quark distribution function q(x) and the fragmentation function D h/q (z). For these "C F -terms" we will here follow the treatment of [22]. The terms with a coefficient N c /2 (we will denote these as the "N c -terms" in the following), i.e. J and J v , have finite transverse integrals 1 but are finite in the limit ξ → 1. Thus they produce a rapidity divergence due to the explicit factor 1/(1 − ξ) in the expression for the multiplicity. If one takes the large N c -limit, there are additional cancellations between some C F and N c -terms that are used in [22]; we will, however, not take this limit here.

A. Explicit subtraction scale
In the CXY calculation the rapidity divergence is subtracted by defining a renormalized dipole cross section as Expressed in coordinate space this reduces to a more familiar looking form in terms of an integral form of the BK renormalization group equation [8,9] S( As pointed out in [27], this subtraction really should include an explicit rapidity factorization scale. The remaining NLO cross section would then depend on this factorization scale, which cancels at this order in perturbation theory against the rapidity up to which the dipole cross section in the leading order cross section is evolved. This is completely analogous with the subtraction of the collinear divergences into DGLAP evolution in the C Fterms, which leaves the NLO cross sections explicitly dependent on a transverse momentum factorization scale.
Let us now accordingly include the rapidity factorization scale and subtract the rapidity divergence as The original choice of [22] simply corresponds to ξ f = 0 in our notation. With this subtraction, and the factorization of the transverse divergence into DGLAP evolution, the cross section now becomes dN pA→hX To see this for the virtual term J v , one has to note that at large q one can use the normalization (9) to perform the l integral, leading to a cancellation between the UV-divergences in the two where P qq is the quark-quark splitting function and the functions I 21 (k ⊥ , ξ) and K(ξ) are defined according to In these expressions we used the plus prescription where f (x) is singular at x = 1 and g(x) is a regular function. Now we can see that if ξ f is chosen to be very close to one (i.e. if one makes sure to subtract only terms that have a very large energy logarithm), the contribution of the last N c -term (integrated from ξ f to 1) is negligible ∼ (1 − ξ f ). The first two N c -terms, on the other hand, yield a large logarithmic contribution ∼ ln(1 − ξ f ) from the upper limit of the integration.

B. Analytical considerations
It is instructive to see how these expressions for the subtracted cross section behave in the opposite limits of small and large transverse momentum for the produced quark. At large transverse momentum the result can be easily obtained from [25], where the unsubtracted cross sections have been matched to collinear perturbation theory. In the large k ⊥ limit only the radiative corrections contribute and the leading behavior comes from where the integrated gluon distribution is with S ⊥ the transverse area of the target hadron.
It is easy to see that if we replace the expression of I in Eq. (1) by its large k ⊥ limit (19), the C F -term will yield a positive contribution. On the other hand, the large k ⊥ behavior of J (20) means that the subtracted N c -term in Eq. (13), will yield a contribution that behaves as a power in k ⊥ and is negative due to the growing ξ-dependence of K in this limit. For values of ξ f close to 0, the magnitude of this negative term is larger than the one of the (positive) C F -term. Therefore, if the leading order cross section falls rapidly at large k ⊥ (as in the GBW parametrization), the leading large k ⊥ behavior comes from the NLO corrections and the cross section is negative. Thus the negativity of the NLO cross section can be traced back to the fact that the unsubtracted cross section at large k ⊥ is proportional to ξ. Subtracting, as in CXY (see (10)), the integral over the whole ξ-interval multiplied by the cross section at ξ = 1 subtracts a very large finite contribution in addition to the divergent part. This oversubtraction of the rapidity divergence is what makes the cross section negative at large transverse momenta. This was the main point in [27], where the authors showed that reintroducing a part of this oversubtracted contribution again makes the cross section positive. When we increase ξ f towards 1, the large negative contribution dies away.
For practical purposes, this is not a formalism that we would wish to use for extremely small transverse momenta, since the independent vacuum fragmentation picture of hadron production is probably not a valid physical picture there. However, from a formal point of view it can be instructive to see what happens in the small k ⊥ limit. The fact that BK evolution preserves the dilute limit S(r) → 1 when r → 0, i.e. the sum rule (9), tells us that the quantity (J (k ⊥ , 1)−J v (k ⊥ , 1)), which is positive in the large k ⊥ limit, must be negative in some other regions of phase space. Indeed, if we denote by δS(k ⊥ ) the contribution that is added to S(k ⊥ ) in one rapidity step of BK evolution, the sum rule ensures that the k-integral is positive. Therefore, if we increase ξ f we subtract a smaller but positive contribution and thus the cross section increases. On the other hand, at smaller k ⊥ we expect (J (k ⊥ , 1) − J v (k ⊥ , 1)) to be negative and the cross section to decrease with increasing ξ f .
As both J and J v are free of IR and UV divergences, we can estimate their leading order behavior by using with |q 0 | > k ⊥ , |q 0 | > l ⊥ , which allows us to figure out the small k ⊥ limit of J and J v as Here Q and Q are some hard momentum scales which are much larger than k ⊥ and q ⊥0 is chosen such that the the above expressions, we see that both J and J v have logarithmic divergences at ξ = 1 and k ⊥ = 0. It is known that (J (k ⊥ , ξ) − J v (k ⊥ , ξ)) should be finite at ξ = 1, and thus the ln(1 − ξ) 2 terms in Eqs. (24) and (25) must have the same coefficient. Therefore S(q ⊥0 ) has to be S(k ⊥ ) at ξ = 1. This also tells us that, at ξ = 1, ln k 2 ⊥ cancels out between the real and virtual terms. Clearly, there is a ln k 2 ⊥ behavior in (J (k ⊥ , ξ) − J v (k ⊥ , ξ)) when ξ = 1. In the GBW model it is possible to study analytically the behavior of the subtraction term. In this model the dipole cross section is given by which leads to For simplicity we will consider here that the saturation scale Q s is a constant. In this model the real N c term reads with where Ei(x) is the exponential integral function, Ei(x) = − ∞ −x dt e −t t , and the virtual term reads Thus the subtraction term is With this expression one can explicitly show that as required to satisfy the sum rule (9). At small k ⊥ , J and J v read and the subtraction term behaves like In Fig. 1, we show the behavior of J (k ⊥ , 1)−J v (k ⊥ , 1) as a function of k ⊥ for a fixed saturation scale Q s = 1 GeV. We observe that, as expected from the general limit (20), the subtraction term is positive at large k ⊥ . The precise  functional form ∼ k 2 ⊥ of the small k ⊥ behavior, on the other hand, is specific to the GBW model. This confirms explicitly that when we increase ξ f , we are subtracting less of the positive contribution in the large k ⊥ region and less of the negative contribution at small k ⊥ . Therefore we would expect that the subtracted multiplicity will increase at large k ⊥ and decrease at small k ⊥ in the GBW model with the increasing of ξ f .

C. Numerical results in the GBW model
In this section we illustrate the features we have discussed so far using the GBW model to parametrize the dipole cross section S(r). Again we focus on the contribution of the quark channel q → q. In this model, the saturation scale Q s at a given x is parametrized by where A is the atomic number of the target and the values of the other parameters are c = 0.56, Q s0 = 1 GeV, x 0 = 3.04 × 10 −4 , λ = 0.288 [30]. In the following the x value at which we evaluate Q s is given by where x g corresponds to the leading order kinematics: Since in the GBW model the dipole cross section has a simple gaussian form, it is possible to perform some integrals analytically. This avoids the need to deal with several oscillatory integrals in the numerical implementation which is thus simpler. It also enhances the difference in the behaviors of the leading and nextto-leading contributions at large k T . Indeed, the leading order contribution is proportional to the Fourier transform of S(r), which is also gaussian, while at large k T the NLO corrections have a power law behavior, as can be seen from Eqs. (19) and (20).
The expression for dσ pA→hX /d 2 pdy h in the GBW model can be found in Ref. [22] in the large N c limit.
Here we need to keep N c finite because we want to have a clear separation between the C F -terms, which are associated with the collinear divergence, and the N c -terms which are associated with the rapidity divergence. In this case, the multiplicity reads dN pA→hX and the expressions for J and J v in the GBW model can be read from Eqs. (28) and (30), respectively. Besides the absence of the impact parameter integration, there are two differences with the corresponding expressions given in Ref. [22]: first, we keep terms proportional to C F − N c /2 which vanish in the large N c limit taken in Ref. [22]. Second, we modify the rapidity divergence subtraction by using the cutoff ξ f introduced previously. In Fig.  2 we show the multiplicity dN pAu→h − X /d 2 pdy h as a function of p ⊥ in the GBW model for ξ f = 0 which corresponds to the choice made in Ref. [22]. We take √ s = 200 GeV, α s = 0.2, µ 2 = 10 GeV 2 and y h = 3.2. For the collinear PDFs q(x) and fragmentation functions D h/q (z) we use the MSTW 2008 [32] and DSS [33] parametrizations, respectively, both at next-to-leading order. As observed in Ref. [26], when ξ f = 0 the NLO multiplicity is negative when p ⊥ is larger than a certain value, of the order of the saturation scale. As discussed above, this negativity comes from the N c -terms. This can be seen from the same figure where we also show the effect of including only the NLO corrections proportional to C F or N c . At large transverse momentum the C F corrections are positive while the N c corrections are negative and large enough to make the total NLO multiplicity negative. We observe that when including only the NLO contributions proportional to N c the multiplicity becomes negative very close to the point where p ⊥ ≈ Q s . To see more clearly the behavior of the NLO corrections at small transverse momentum, we show in Fig. 3 the ratio of the NLO and LO multiplicity for p ⊥ ≤ 2.5 GeV when including only the C F or N c contributions or both. One can note that, already for values of p ⊥ of the order of 2 GeV, both the C F and N c contributions are of the same order of magnitude as the LO term. As discussed in Sec. III, instead of subtracting the whole ξ interval when subtracting the rapidity divergence (10), one can introduce an explicit cutoff ξ f to  determine which contributions are included in the renormalized dipole cross section (12). We recall that, compared to the results shown above for ξ f = 0, only the N c -terms are affected by this procedure. In Fig. 4, we show the multiplicity for several fixed values of ξ f as a function of p ⊥ . We observe that, for the reason exposed in Sec. III B, values of ξ f close to 1 lead to a positive multiplicity up to larger values of p ⊥ . In particular, for ξ f 0.999, the multiplicity is positive up to p ⊥ = 8 GeV. In Fig. 5, we show the ratio NLO/LO for several values of ξ f and p ⊥ ≤ 2.5 GeV. Here we see that values of ξ f very close to 1, corresponding to a positive multiplicity at large p ⊥ , lead to a NLO multiplicity smaller than at leading order at moderate p ⊥ . In Fig. 6, we show the same ratio with a fixed value of the saturation scale, Q s = 1 GeV. Here we see clearly that, as could be expected from the small k ⊥ behavior of the subtraction term (35), larger values of ξ f lead to smaller cross sections at small p ⊥ and larger cross sections at large p ⊥ . To conclude this section, it appears that the choice of the value of ξ f can have an important impact on the final results, both at small and large transverse momentum. In the following we propose a way to fix the value of this parameter based on physical considerations.

IV. ORDERING IN LIGHT CONE ENERGY
It has recently been emphasized [31] that, in order to have a stable BK evolution beyond leading order, the gluon cascade resummed by the evolution should be ordered in the light cone energy k − of the projectile. The requirement of ordering in k − can equivalently be thought of as an ordering in its conjugate variable x + , which is the light cone lifetime or "Ioffe time" of fluctuations in the projectile [28]. In many works, this feature is known as Figure 7. Gluon emission the "kinematical constraint" or as imposing "exact kinematics" [25,29,34], a terminology that we will comment on in Sec.V. Imposing k − -ordering has also been one ingredient in the program of "small-x resummation" in the context of the linear BFKL evolution [35][36][37][38]. Since we are working in a frame where we consider the gluons as being emitted from the probe, the emitted gluon always has a smaller k + momentum than its parent quark; the emissions are therefore naturally ordered in k + . The requirement of k − ordering, on the other hand, must be imposed separately. At leading order, the physical picture of the scattering is that of an incoming collinear quark (i.e. with only a k + momentum component) that acquires a transverse momentum and light cone energy (k and k − ) from the target. Considering the production of a quark with a fixed k, the light cone energy required is The fraction of the target momentum required to set the produced quark on shell is defined by the leading order kinematics as x g = k − LO /P − , where P − is the momentum of the target. Thus we can write 1/(2x p P + ) = x g P − /k 2 .
We now want to implement k − -ordering in the calculation of the single inclusive cross section. For this purpose let us consider the emission of a gluon from the incoming quark that is present in all of the contributing diagrams. Labeling the transverse momentum of the incoming quark as q and radiated gluon as l (cf. Fig. 7), the light cone energy introduced from the gluon emission is where 1 − ξ is the momentum fraction of the emitted gluon. In some diagrams, the splitting happens before the interaction with the target and q = 0, in final state radiation diagrams q ⊥ ∼ Q s . Similarly, for some diagrams, but not all, the quark does not interact with the target after the emission and q − l = k. In general, the momenta l, q in Eq. (40) are integrated over in the diagram, whereas the transverse momentum k is that of the produced quark and is kept fixed (at a fixed fragmentation z). It appears as an overall prefactor in (40) when the incoming quark longitudinal momentum x p P + is expressed in terms of the target variables x g and P − and is therefore common to all diagrams contributing to the cross section. We now want to derive a renormalization group equation describing the target. This means that emissions in a certain kinematical regime have to be subtracted from the cross section and absorbed into a redefinition of the target, as was done in Eq. (10). In particular, this kinematical regime to be subtracted should include the limit ξ = 1. Now the requirement of a k − -ordering in the evolution means that the subtraction criterion should be a condition on ∆k − . Since ∆k − ∼ 1/(1 − ξ) for ξ → 1, this means that fluctuations around ξ = 1 with ∆k − larger than a certain factorization scale should be subtracted. Thus, our renormalization scheme is defined by the condition that all the contributions with ξ ≥ ξ f (k ⊥ ) are to be subtracted from the cross section. Most importantly, we want the subtracted region to include contributions with where x f is in principle an arbitrary factorization scale that will appear both in the leading order term as the rapidity up to which the target must be evolved, and in the hard factors of the NLO cross section. We emphasize that (41) is not a (momentum conservation) kinematical constraint for the process, but a renormalization condition specifying which parts of phase space should be subtracted from the cross section. A natural choice that resums all the large energy logarithms into the evolution is to take x f ≈ x g . In the typical kinematical regime when all the transverse momenta are of the same order, ∆k − qg ∼ x g /(1 − ξ) x f for all ξ, and we can safely subtract the whole ξ interval as done in [21,22]. However, when the produced quark momentum k is much larger than the typical target scale Q s , there are configurations where the k − -ordering condition (41) is not satisfied for all ξ, because l is integrated over in a range that includes typical target scales l ⊥ ∼ Q s . Thus we should, for large k, only subtract values of ξ close to 1 that satisfy This line of argument leads us to our proposal to implement k − -ordering or the "kinematical constraint" by subtracting a k ⊥ -dependent fraction of the ξ integral from the cross section as Note that since there is a factorization scale x g /x f in the limit, we have the freedom to change this scale with a factor of order 1. This change should cancel against a corresponding change in the rapidity to which the dipole amplitude is evolved. An alternate, somewhat smoother form with the same parametric behavior would be where we can vary the factorization scale x g /x f in an interval such as 1 2 . . . 2. In practice, the two choices (45) and (46) lead to very similar results. Therefore we will only show results obtained with the smoother choice (46). In Figs. 8 and 9 we show the multiplicity as a function of p ⊥ for three values of the ratio x g /x f . In the three cases, the NLO multiplicity still becomes negative for some value of the transverse momentum. For x g /x f = 1 or 2 this happens approximately at the same point as for ξ f = 0, as can be seen by comparison with Fig. 2. On the other hand, for x g /x f = 0.5 the multiplicity is positive up to a significantly larger value of p ⊥ of the order of 6 GeV. For comparison we also show in Fig. 8 the results obtained by following the approach of Ref. [29], where the original CXY subtraction is used but a kinematical constraint is imposed which leads to an additional correction term L q . This additional term extends the p ⊥ range where the multiplicity is positive by about 0.5 GeV, which is similar to what was obtained in Ref. [29] in the same kinematics. This is the main result of our paper: we have shown that it is possible to choose a value of x g /x f which is still in its "natural" range and extends significantly the region where the NLO cross section has a reasonable physical behavior. In the case we have studied here the dependence on the exact choice for the renormalization scale is stronger than one would like, but we believe that this is due to two aspects of the calculation that can be improved in the future. First, we tried here to implement ordering in an effective way by relying only on external scales to fix ξ f . For more accurate results, one should instead impose the ordering (41) inside the transverse momentum integrals in J and J v . This could lead to sizeable differences with the results we have shown since, as can be seen from Fig. 8, our results are still very sensitive to the value of ξ f . Second, we have used the GBW model to parametrize the dipole cross section. While this has practical advantages like enabling us to perform some integrals analytically, in this model the leading and nextto-leading order contributions have very different behaviors at large transverse momentum: the leading order term behaves like a gaussian, while the next-to-leading order corrections behave like a power law. Therefore, at x f x g compared with the results at leading order and when using the correction term L q introduced in Ref. [29].  Figure 9. Ratio of the multiplicity at next-to-leading and leading order using the parametrization (46) of ξ f for different values of large transverse momentum, the behavior of the multiplicity is governed entirely by the NLO corrections, which is quite unnatural. On the other hand, using a more physical dipole cross section, such as one obtained by solving the Balitsky-Kovchegov [8,9] equation, should make the leading order contribution behave more like a power law. This means that the end result would be less sensitive to the exact choice made for the parametrization of ξ f . Taking a BK-evolved dipole cross section instead of the GBW one would also make the cross section formally independent of the factorization scale x f at this order in α s , which should significantly reduce the factorization scale dependence.

V. DISCUSSION AND OUTLOOK
To summarize, we have in this paper proposed to modify the subtraction procedure of the rapidity divergence in the calculation of single inclusive hadron production in the hybrid formalism [21,22]. It has been hoped that imposing ordering in the light cone energy of the probe could solve the issue of large negative NLO corrections to the cross section. We have here argued that the most natural way to implement this ordering is to impose it on the kinematics of the part that is subtracted from the cross section and absorbed into the target. To explicitly demonstrate the effect of our proposal, we have performed the calculations at finite N c for the GBW model. The part corresponding to the rapidity divergence is proportional to N c , whereas the parts associated with DGLAP evolution have a color factor C F . Compared to the CXY formulation, our suggestion only modifies the N c -terms and leaves the C F -terms unchanged. Indeed, we have shown that such a modification of the subtraction scale can lead to a more stable perturbative expansion for this cross section at high p ⊥ . By a suitable choice of the factorization scale x f , the stability of the perturbative expansion can be extended to arbitrarily high p ⊥ . This should be contrasted with the recent calculation of Ref. [29], where a similar correction is obtained as an additional correction which in some cases still leaves the cross section negative at high enough p ⊥ . The main difference between our approach here and that of Ref. [29] is that there the kinematical modification is treated as a single correction term whereas here it is resummed to all orders by shifting the evolution variable in the BK equation.
In some works (e.g. [25,29]), the kinematical constraint is explained as a requirement of keeping the momentum fraction in the target x a as less than 1. We would rather prefer to characterize this requirement by saying that the x a of the corrections that are resummed into the BK evolution should be larger than the x g determined by the leading order kinematics. The whole formalism here is based on the eikonal approximation, which is valid for only small enough x a . If the actual result for the cross section were to really depend on the dipole cross section at very large target x a , the whole formalism would be in grave trouble. Fortunately this is not the case, since the contribution from large x a ∼ 1 (corresponding to (1 − ξ) x g ) is subtracted from the cross section and absorbed into the renormalization group evolution of the target: in fact it is only a tiny part of the subtraction. In spite of this different terminology, the actual equations in Refs. [25,29] lead to the same parametric dependence that we have here: the longitudinal factorization scale has to be modified by the presence of a large transverse momentum logarithm ∼ ln k 2 ⊥ /Q 2 s . In order to perform the whole calculation consistently at NLO accuracy, one needs to also solve the NLO version of the BK equation. To do this, the actual BK equation to be solved must be consistent with the subtrac-tion procedure in the cross section. Different subtraction schemes correspond to different versions of the BK equation. Since the rapidity dependence of the dipole cross section is proportional to α s , a modification of the subtraction scale in Eq. (44) is a higher order α 2 s effect, and the difference appears at the NLO level for the BK equation. In fact, it was recently shown in [16] that a problematic double logarithmic correction to the NLO BK equation [11] can indeed be resummed by redefining the evolution variable in a similar way as we have done in Eq. (44). When written as a differential equation in Y = ln 1/x f , the modified BK equation is nonlocal in rapidity, which is inconvenient for a numerical solution. The authors of [16] proposed a clever way of rewriting this in a local form, which has successfully been solved numerically [17,18].
Compared to the work in [16][17][18], our proposal here has been a very simplistic one. By explicitly taking Q s as the scale of the logarithm, we have been able to write a subtraction procedure in such a way that we can reuse most of the CXY formulae with only slight modifications. To be consistent [16][17][18], the momentum scale in Eq. (45) should actually be one that is integrated over in J and J v . This would require a new computation of the regularized "hard factors" remaining after the subtraction. For a calculation that is consistent with recent NLO BK solutions it would be useful to carry out this subtraction for the resummation proposed in [16]. This is, however, left for future work.