On the implementation of NLO high energy factorization in single inclusive forward hadron production

TL;DR

This paper investigates the implementation of NLO high energy factorization in single inclusive forward hadron production, demonstrating stability and positivity of the cross section with fixed and running coupling, and addressing issues in previous calculations.

Contribution

It introduces a new rapidity factorization procedure for NLO calculations, ensuring physically meaningful and stable cross sections in forward hadron production.

Findings

Fixed coupling yields positive, stable cross sections.

Running coupling results are stable but differ from previous coordinate space approaches.

The approach improves theoretical consistency in high energy QCD calculations.

Abstract

Single inclusive particle production cross sections in high energy hadron collisions at forward rapidity are an important benchmark process for the CGC picture of small x QCD. Recent calculations of this process have not led to a stable perturbative expansion for this quantity at high transverse momenta. We consider the quark channel production cross section using the new rapidity factorization procedure proposed by Iancu et al. We show that for fixed coupling one does indeed obtain a physically meaningful cross section which is positive and reduces in a controlled way to previous leading order calculations. We also consider a running coupling that depends on the transverse momentum of the produced particle. This gives a stable result which, however, is not fully consistent with previous leading order calculations that use a coordinate space running coupling.

Figures

Equations

$$\mathcal{S}(k_{\perp })=\mathcal{S}(k_{\perp },\mathbf{b})=\int {\mathrm{d}}^2\mathbf{r}{e}^{-i\mathbf{k}\cdot \mathbf{r}}S(\mathbf{r}),$$

$$S(\mathbf{r}=\mathbf{x}-\mathbf{y})=⟨\frac{1}{N_{\mathrm{c}}}\mathrm{Tr}V(\mathbf{x})V^{†}(\mathbf{y})⟩,$$

$$\frac{\mathrm{d}{N}^{pA\to qX}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}=x_{p}q(x_p)\frac{{\mathcal{S}}_0(k_{\perp })}{(2\pi {)}^2}$$

$$-\frac{{\alpha }_{\mathrm{s}}}{2{\pi }^2}{\int }_0^{{\xi }_{\text{max}}}\mathrm{d}\xi \frac{1+{\xi }^2}{1-\xi }{x}_{p}q\left(x_p\right)\left\{C_{\mathrm{F}}{\mathcal{I}}_v(k_{\perp },\xi ,X(\xi ))+\frac{N_{\mathrm{c}}}{2}{\mathcal{J}}_v(k_{\perp },\xi ,X(\xi ))\right\}.$$

$$\mathcal{I}(k_{\perp },\xi ,X(\xi ))$$

$$\mathcal{J}(k_{\perp },\xi ,X(\xi ))$$

$${\mathcal{I}}_v(k_{\perp },\xi ,X(\xi ))$$

$${\mathcal{J}}_v(k_{\perp },\xi ,X(\xi ))$$

$$\frac{\mathrm{d}{N}^{\text{LO}+N_{\mathrm{c}}}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}=$$

$$+{\alpha }_{\mathrm{s}}{\int }_0^{1-x_g/x_0}\frac{\mathrm{d}\xi }{1-\xi }\mathcal{K}(k_{\perp },\xi ,X(\xi )),$$

$$\mathcal{K}(k_{\perp },\xi ,X(\xi ))=$$

$$\displaystyle\times\bigg{[}$$

$$\displaystyle-x_{p}q\left(x_{p}\right)\mathcal{J}_{v}(k_{\perp},\xi,X(\xi))\bigg{]}.$$

$$\frac{\mathrm{d}{N}^{\text{LO}+N_{\mathrm{c}}}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}=$$

$$+{\alpha }_{\mathrm{s}}{\int }_0^{1-x_g/x_0}\frac{\mathrm{d}\xi }{1-\xi }\mathcal{K}(k_{\perp },\xi ,X(\xi ))$$

$$\equiv$$

$$\mathcal{S}(k_{\perp },x_g)=\mathcal{S}(k_{\perp },x_0)$$

$$+2{\alpha }_{\mathrm{s}}{N}_{\mathrm{c}}\underset{0}{\overset{1-x_g/x_0}{\int }}\frac{\mathrm{d}\xi }{1-\xi }\left[\mathcal{J}(k_{\perp },1,X(\xi ))-{\mathcal{J}}_v(k_{\perp },1,X(\xi ))\right],$$

$$\begin{gathered} \frac{\mathrm{d}{N}^{\text{LO}+N_{\mathrm{c}}}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}=x_{p}q(x_p)\frac{\mathcal{S}(k_{\perp },x_g)}{(2\pi {)}^2} \\ +{\alpha }_{\mathrm{s}}\underset{0}{\overset{1-x_g/x_0}{\int }}\frac{\mathrm{d}\xi }{1-\xi }\left[\mathcal{K}(k_{\perp },\xi ,X(\xi ))-\mathcal{K}(k_{\perp },1,X(\xi ))\right] \\ \equiv \frac{\mathrm{d}{N}^{\text{LO}}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}+\frac{\mathrm{d}{N}^{N_{\mathrm{c}},sub}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}. \end{gathered}$$

$$\begin{gathered} \frac{\mathrm{d}{N}^{\text{LO}+N_{\mathrm{c}}}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}=x_{p}q(x_p)\frac{\mathcal{S}(k_{\perp },x_g)}{(2\pi {)}^2} \\ +{\alpha }_{\mathrm{s}}{\int }_0^{1-x_g/x_0}\frac{\mathrm{d}\xi }{1-\xi }\left[\mathcal{K}(k_{\perp },\xi ,x_g)-\mathcal{K}(k_{\perp },1,x_g)\right], \end{gathered}$$

$$\frac{\mathrm{d}{N}^{C_{\mathrm{F}}}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}\equiv \frac{{\alpha }_{\mathrm{s}}}{2{\pi }^2}{C}_{\mathrm{F}}\left[\underset{x_p}{\overset{1-x_g/x_0}{\int }}\mathrm{d}\xi \frac{1+{\xi }^2}{1-\xi }\frac{x_p}{\xi }q\left(\frac{x_p}{\xi }\right){\mathcal{I}}^{\text{finite}}(k_{\perp },\xi ,X(\xi ))-\underset{0}{\overset{1-x_g/x_0}{\int }}\mathrm{d}\xi \frac{1+{\xi }^2}{1-\xi }{x}_{p}q\left(x_p\right){\mathcal{I}}_v^{\text{finite}}(k_{\perp },\xi ,X(\xi ))\right],$$

$${\mathcal{I}}^{\text{finite}}(k_{\perp },\xi ,X(\xi ))$$

$${\mathcal{I}}_v^{\text{finite}}(k_{\perp },\xi ,X(\xi ))$$

$$S(\mathbf{r},x_0)=\exp \left[-\frac{{\mathbf{r}}^2{Q}_{\mathrm{s},0}^2}{4}\ln \left(\frac{1}{|\mathbf{r}|{\Lambda }_{\mathrm{QCD}}}+e\right)\right],$$

$$\frac{\mathrm{d}{N}^{\text{NLO}}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}=\frac{\mathrm{d}{N}^{\text{IC}}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}+\frac{\mathrm{d}{N}^{N_{\mathrm{c}},unsub}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}+\frac{\mathrm{d}{N}^{C_{\mathrm{F}}}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}.$$

$$\frac{\mathrm{d}{N}^{\text{NLO}}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}=\frac{\mathrm{d}{N}^{\text{LO}}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}+\frac{\mathrm{d}{N}^{N_{\mathrm{c}},sub}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}+\frac{\mathrm{d}{N}^{C_{\mathrm{F}}}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}.$$

$$\frac{\mathrm{d}{N}^{\text{LO}+C_{\mathrm{F}}}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}=\frac{\mathrm{d}{N}^{\text{LO}}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}+\frac{\mathrm{d}{N}^{C_{\mathrm{F}}}}{{\mathrm{d}}^2\mathbf{k}\mathrm{d}y}.$$

$$S(\mathbf{r},x_0)=\exp \left[-\frac{{\mathbf{r}}^2{Q}_{\mathrm{s},0}^2}{4}\ln \left(\frac{1}{|\mathbf{r}|{\Lambda }_{\mathrm{QCD}}}+e_c\cdot e\right)\right],$$

$${\alpha }_{\mathrm{s}}({\mathbf{r}}^2)=\frac{4\pi }{{\beta }_0\ln \left(\frac{4C^2}{{\mathbf{r}}^2{\Lambda }_{\mathrm{QCD}}^2}\right)},$$

$${\alpha }_{\mathrm{s}}(k_{\perp }^2)=\frac{4\pi }{{\beta }_0\ln \left(\frac{C_{\text{mom}}^2{k}_{\perp }^2}{{\Lambda }_{\mathrm{QCD}}^2}\right)}$$

$${\mathcal{I}}^{\text{finite}}(k_{\perp },\xi ,X(\xi ))=\int {\mathrm{d}}^2\mathbf{r}$$

$$\displaystyle-2e^{-i{\mathbf{k}}\cdot{\mathbf{r}}}\int\frac{\mathrm{d}^{2}{\mathbf{x}}}{(2\pi)^{2}}\frac{{\mathbf{x}}\cdot({\mathbf{x}}+{\mathbf{r}})}{{\mathbf{x}}^{2}({\mathbf{x}}+{\mathbf{r}})^{2}}S(\xi{\mathbf{r}}-(1-\xi){\mathbf{x}},X(\xi))\bigg{]}\,,$$

$${\mathcal{I}}_v^{\text{finite}}(k_{\perp },\xi ,X(\xi ))=\int {\mathrm{d}}^2\mathbf{r}$$

$$\mathcal{J}(k_{\perp },\xi ,X(\xi ))=\int {\mathrm{d}}^2\mathbf{r}$$

$$=\int {\mathrm{d}}^2\mathbf{r}$$

$${\mathcal{J}}_v(k_{\perp },\xi ,X(\xi ))=\int {\mathrm{d}}^2$$

$$=\int {\mathrm{d}}^2\mathbf{r}$$

$$\partial_{Y}S({\mathbf{r}},X)=2\alpha_{\mathrm{s}}{N_{\mathrm{c}}}\left[\widetilde{\mathcal{J}}({\mathbf{r}},1,X)-\widetilde{\mathcal{J}}_{v}({\mathbf{r}},1,X)\right]=2\alpha_{\mathrm{s}}{N_{\mathrm{c}}}\int\frac{\mathrm{d}^{2}{\mathbf{x}}}{(2\pi)^{2}}\frac{{\mathbf{r}}^{2}}{{\mathbf{x}}^{2}({\mathbf{x}}+{\mathbf{r}})^{2}}\big{[}S(-{\mathbf{x}},X)S({\mathbf{r}}+{\mathbf{x}},X)-S({\mathbf{r}},X)\big{]},$$

$${\partial }_{Y}S(\mathbf{r},X)=2{\alpha }_s({\mathbf{r}}^2)N_{\mathrm{c}}\int \frac{{\mathrm{d}}^2\mathbf{x}}{(2\pi {)}^2}$$

Full text

I Introduction

To study the nonlinear gluon saturation regime one needs to look at processes that involve a momentum fraction $x$ in the target as small as possible, while still having a large enough transverse momentum to justify a weak coupling treatment. At LHC energies a good process to achieve this is particle production at forward rapidities. In the CGC picture this process can be calculated using the “hybrid” formalism, where a quark or a gluon, taken from the usual collinear parton distribution in the probe at large $x$, passes through the target color field. It receives a kick from the intrinsic transverse momentum of the target gluons, leading to the $p_{\perp}$ spectrum of the produced hadrons.

Recently much work has been done to calculate forward particle production in the hybrid picture to next-to-leading order (NLO) accuracy in the QCD coupling. Cross sections for this process were calculated to NLO accuracy in Refs. Chirilli:2011km ; Chirilli:2012jd (see also the earlier works Dumitru:2005gt ; Altinoluk:2011qy ; JalilianMarian:2011dt ), where soft and collinear divergences in the one-loop calculation were factorized respectively into the BK Balitsky:1995ub ; Kovchegov:1999yj evolution of the target and the DGLAP evolution of the probe proton and the fragmentation functions. In the first numerical implementation Stasto:2013cha of the “CXY” factorization framework of Chirilli:2011km ; Chirilli:2012jd the NLO corrections turned out to be large and render the total cross section negative at large transverse momenta of the produced particle, signaling a problem in the organization of the perturbative series. Following this observation, several interpretations of the origin of the problem have been proposed Kang:2014lha ; Altinoluk:2014eka ; Watanabe:2015tja .

We recently argued Ducloue:2016shw that the origin of the problem lies in the subtraction procedure where the soft divergence in the NLO calculation is factorized into the BK evolution of the target. Following this discussion, Iancu et al. Iancu:2016vyg suggested a new formulation of the NLO cross section in a way that explicitly yields a positive cross section. The formulation is based on the observation that in a certain limit which we will discuss in more detail below, the expression for the cross section involves similar structures as an integral form of the BK equation. If one takes care not to break this equivalence by further approximations and chooses a Fourier-positive initial condition for the BK equation, this will guarantee the positivity of the cross section. The resulting expression can either be written as a “subtracted” cross section, or in a form where no explicit subtraction is performed.

The purpose of this paper is to present a practical numerical implementation of the manifestly positive formulation for the cross section presented in Iancu:2016vyg . We will first briefly review the expressions for the cross section in Sec. II. We then in Sec. III explicitly show that this procedure works at fixed coupling as anticipated. In Sec. IV we discuss the problems associated with introducing a running coupling into the expression. We implement the calculation using a running coupling that depends on the transverse momentum of the produced particle, which would be the natural thing to do for a cross section. We show that since the BK equation, and fits to DIS data using it, are usually implemented in coordinate space, this introduces a mismatch between the manifestly positive form of the cross section and a “subtracted” version where the separation between LO and NLO contributions is explicit. We then conclude in Sec. V with a brief outlook for practical phenomenological applications of the formalism. In Appendix A we discuss an alternative coordinate space formulation of the cross section, which does not suffer from this mismatch between manifestly positive and subtracted formulations, but results in unphysically large values for the NLO corrections to the cross section.

II Single inclusive particle production at NLO

Our starting point are the CXY formulae derived in Refs. Chirilli:2011km ; Chirilli:2012jd . We will concentrate here on the quark channel, for which the leading order cross section is proportional to the Fourier-transform of the dipole operator

[TABLE]

where $V({\mathbf{x}})$ is a fundamental representation Wilson line in the color field of the target. The NLO cross section involves also adjoint representation Wilson lines from the gluon interacting with the target color field. Using Fierz identities and the mean field approximation, which replaces expectation value of products of dipole operators by products of expectation values, also the NLO cross section can be expressed in terms of $\mathcal{S}(k_{\perp})$. Since our goal is to study the negativity problem and how it would be affected by the proposal of Ref. Iancu:2016vyg , we leave out the fragmentation functions which do not play any role here. Following the notation in Ducloue:2016shw , the (unsubtracted) CXY quark multiplicity can be written as:

[TABLE]

Here we have left $\xi_{\text{max}}$, the upper limit for the $\xi$-integrals, unspecified for now as we will return to this question shortly. The kinematical variables are defined as $x_{p}=k_{\perp}e^{y}/\sqrt{s}$, $x_{g}=k_{\perp}e^{-y}/\sqrt{s}$ and $k_{\perp}=|{\mathbf{k}}|$. The most important one for our discussion here is the momentum fraction $\xi$: the fragmenting quark carries a fraction $\xi$ of the incoming quark longitudinal momentum. Thus the incoming quark has a momentum fraction $x_{p}/\xi$ of the incoming proton, where $x_{p}$ is the probe momentum fraction in the leading order kinematics. The radiated gluon in the NLO terms carries a longitudinal momentum fraction $1-\xi$, i.e. the limit $\xi\to 1$ corresponds to the soft gluon emission that must be resummed into the BK evolution of the target. For the following discussion it is important to note the interpretation of the variable $x_{g}$: it is easy to see that for producing a final state quark with transverse momentum $k_{\perp}$ at leading order, $x_{g}$ is the fraction of target longitudinal momentum $P^{-}$ needed to put this quark (with longitudinal momentum $x_{p}P^{+}$) on shell. Thus, in a leading order calculation such as Lappi:2013zma , with BK evolution starting from an initial momentum fraction $x_{0}\sim 0.01$, one would evolve the target for $\ln x_{0}/x_{g}$ units in rapidity (here rapidity is defined as the logarithm of the inverse of the target momentum fraction).

The dependence on the Wilson line correlators in the target comes through transverse momentum integrals that can be expressed in terms of the dipole operator $\mathcal{S}$:

[TABLE]

At this stage we have not specified how the dipole operators depend on the energy or rapidity, but have just denoted this dependence by $X(\xi)$. This dependence is indeed the crux of the argument in Ref. Iancu:2016vyg , which we will here reproduce in a slightly altered form. As emphasized e.g. in Altinoluk:2014eka ; Ducloue:2016shw , one achieves a more stable perturbative expansion when this energy dependence is related to the momentum fraction in the target, i.e. the $k^{-}$ scale in the scattering. This distinction becomes important at large transverse momenta for the produced particle Kang:2014lha ; Altinoluk:2014eka ; Watanabe:2015tja ; Ducloue:2016shw . However, the formulation of Ref. Iancu:2016vyg allows one to get a physically reasonable cross section even without imposing a separate “kinematical constraint,” within just the usual “Regge” kinematics, where all transverse momenta are assumed to be of the same order. Thus we will for now keep this assumption, which allows us to relate the momentum fraction in the target to the kinematics of the probe as $X(\xi)\approx x_{g}/(1-\xi)$.

The CXY cross sections have been derived using the eikonal approximation, which allows one to describe the interaction with the target in terms of Wilson lines. This approximation is only valid for sufficiently high energy scattering, where now “high energy” should refer to the whole quark-gluon state. This is usually reflected in a restriction that the momentum fraction $X(\xi)$ at which one evaluates the dipole cross sections must be smaller than some limiting value $x_{0}\sim 0.01$. For high energy evolution this corresponds to a nonperturbative initial condition for the evolution of the target which is fit to experimental data. With $X(\xi)=x_{g}/(1-\xi)$ this means that the $\xi$-integration must be restricted to $\xi<1-x_{g}/x_{0}\equiv\xi_{\text{max}}$. This now defines the upper limit of the integration that was left unspecified in Eq. (II). Note that this restriction implies that extrapolating beyond the physical region $\xi_{\text{max}}>x_{p}$ to the kinematical limit $x_{g}=x_{0}$, the NLO corrections are explicitly set to approach zero. Calculating NLO corrections involving larger target longitudinal momenta $k^{-}$ is strictly speaking not possible in this formalism and would require going beyond the eikonal approximation. We would argue that setting these contributions to zero is a more controlled approximation than letting them evaluate to some arbitrary value. This should be kept in mind when comparing these calculations to experimental data: this restriction means that when one approaches the kinematical limit $x_{g}\to x_{0}$, the phase space for NLO contributions in our calculation is cut off by this constraint. In this limit the calculation could eventually be matched onto collinear factorization, which becomes more appropriate when $x_{g}$ is large, as done in Ref. Stasto:2014sea .

II.1 ${N_{\mathrm{c}}}$-terms

In Ref. Ducloue:2016shw it was shown that the negativity of the CXY cross section at large transverse momentum is caused by the NLO corrections proportional to ${N_{\mathrm{c}}}$, which become very large and negative while the NLO corrections proportional to $C_{\mathrm{F}}$ are positive. Therefore we first consider only the leading order and ${N_{\mathrm{c}}}$-terms, as was done in Ref. Iancu:2016vyg . These contributions can be written as

[TABLE]

where $\mathcal{K}$ is defined as

[TABLE]

At this stage the dipole operator $\mathcal{S}_{0}$ is a “bare” one and does not evolve with rapidity. Since $\mathcal{K}(k_{\perp},\xi,X)$ approaches a nonzero value for $\xi\to 1$ at fixed $X$, the cross section contains a large logarithmic integral in the high energy limit $x_{g}\to 0$; this should be resummed using the BK equation. The “bare” dipole $\mathcal{S}_{0}$ is then identified with the initial condition of the BK evolution, formulated at a given initial rapidity $\ln\frac{1}{x_{0}}$. We thus rewrite (II.1) as

[TABLE]

where “IC” stands for the multiplicity at the initial rapidity scale $x_{0}$. As long as the initial dipole amplitude is Fourier-positive (which is a nontrivial requirement Giraud:2016lgg that we will however assume to be satisfied here, as it is for the parametrizations that we will use), the multiplicity (II.1) is positive up to large transverse momenta. Now, since we are taking $X(\xi)=x_{g}/(1-\xi)$ as is appropriate in the Regge kinematics, we have $\mathrm{d}\ln X(\xi)=-\mathrm{d}\ln(1-\xi)$. This enables us to write an integral version of the BK equation as

[TABLE]

and to write, assuming that $x_{p}<1-x_{g}/x_{0}$ ,

[TABLE]

Since this equation is equivalent to (II.1), it is positive even up to large transverse momenta. Equation (13) is written explicitly as a sum of a leading order contribution (with a BK-evolved dipole) and a contribution proportional to $\alpha_{\mathrm{s}}$ that has no large energy logarithm because $\mathcal{K}(k_{\perp},\xi,X(\xi))-\mathcal{K}(k_{\perp},1,X(\xi))$ vanishes for $\xi\to 1$. Thus it can be naturally interpreted as the sum of a LO contribution and a NLO correction. In the version (II.1), on the other hand, the way to reduce the NLO expression to the LO limit is less transparent. Instead of dropping out a term explicitly proportional to $\alpha_{\mathrm{s}}$, the LO limit of (II.1) is taken by replacing $\mathcal{K}(k_{\perp},\xi,X(\xi))$ with $\mathcal{K}(k_{\perp},1,X(\xi))$, i.e. setting $\xi=1$ inside the kernel, but without changing the rapidity scale of the dipole correlators.

To go from (II.1) or (13) to the CXY expressions one starts by noting that, because of the subtraction $\mathcal{K}(k_{\perp},\xi,X(\xi))-\mathcal{K}(k_{\perp},1,X(\xi))$, the integral (13) is dominated by the lower limit $\xi\ll 1$. Thus one could argue that the rapidity of the dipole amplitude can be replaced by its value at $\xi=0$, i.e. $\mathcal{S}(k_{\perp},X(\xi))$ by $\mathcal{S}(k_{\perp},X(0))=\mathcal{S}(k_{\perp},x_{g})$ and similarly for $\mathcal{S}(q_{\perp},X(\xi))$ and $\mathcal{S}(l_{\perp},X(\xi))$. This approximation

[TABLE]

is perfectly justified in a weak coupling sense. However, as pointed out in Ducloue:2016shw ; Iancu:2016vyg , it makes the cross section negative, because at large transverse momentum the subtracted NLO term $\mathcal{K}(k_{\perp},\xi,x_{g})-\mathcal{K}(k_{\perp},1,x_{g})$ is negative and can dominate over the leading order term. Thus the necessary ingredient in keeping the result physically meaningful turns out to be to not just evaluate the dipole amplitude at the rapidity scale $x_{g}$ of the leading order cross section, but to keep the dependence of $X(\xi)$ on the integration variable.

In addition to changing the rapidity argument of the dipole operators, another approximation is needed to recover the CXY subtraction result, namely to replace $1-x_{g}/x_{0}$ in the upper limit of the integration over $\xi$ in Eq. (14) by 1. Now that the rapidity argument $X(\xi)$ in the dipole operators has been replaced by $x_{g}$, changing the integration limit is formally possible without extending the dipole parametrization to the large $x$ region where the eikonal approximation is not valid. But this would make the problem disappear only superficially: any contribution from $\xi>1-x_{g}/x_{0}$ would come from the region where the derivation of the original expression (II) is dubious because of the large invariant mass of the produced quark-gluon system. The question of the proper value of the upper limit in $\xi$, or of the correct dependence $X(\xi),$ for large transverse momenta of the produced particle is the “kinematical constraint” or “Ioffe time” issue that has been extensively discussed by several authors Altinoluk:2014eka ; Watanabe:2015tja ; Ducloue:2016shw . We will however not attempt to take these corrections into account here, but stay within the Regge kinematical approximation $X(\xi)=x_{g}/(1-\xi)$.

In the rest of this paper we will refer to the explicitly positive formulation (II.1) as the “unsubtracted” and its subtracted version (13) as the “subtracted” formulation.

II.2 $C_{\mathrm{F}}$-terms

Let us now turn to the $C_{\mathrm{F}}$-terms. While they do not pose similar manifest problems as the ${N_{\mathrm{c}}}$-terms, one nevertheless has to make certain choices concerning the rapidity dependence of the dipole correlators when evaluating the cross section. Writing explicitly the integration limits and subtracting the $1/\varepsilon$ poles corresponding to the DGLAP evolution of the probe quark distributions and fragmentation functions (which are not written explicitly here), the $C_{\mathrm{F}}$-terms become

[TABLE]

where $\mathcal{I}^{\text{finite}}$ and $\mathcal{I}_{v}^{\text{finite}}$, obtained after subtracting the collinear divergence from $\mathcal{I}$ and $\mathcal{I}_{v}$, read

[TABLE]

To arrive at the expressions (15)-(17) we have adopted a specific choice for the previously mentioned ambiguity: namely the rapidity scale at which the dipole operators are evaluated in the $C_{\mathrm{F}}$-terms. To be more precise: to get a finite result for the $C_{\mathrm{F}}$-terms one must subtract from the cross section collinear infrared divergences that manifest themselves as $1/\varepsilon$ poles in the cross section. Thus one is subtracting terms that are proportional to the bare dipole operator $\mathcal{S}_{0}$ from the NLO cross section (II) and neglecting a collinearly divergent remainder $\sim\alpha_{\mathrm{s}}C_{\mathrm{F}}(1/\varepsilon)\left(\mathcal{S}(X(\xi))-\mathcal{S}_{0}\right)$. Such terms are formally of NNLO order, since the difference $\left(\mathcal{S}(X(\xi))-\mathcal{S}(x_{0})\right)$ is proportional to $\alpha_{\mathrm{s}}$. This difference can, however, be numerically large, indeed the large difference between $\mathcal{S}(X(\xi))$, $\mathcal{S}(x_{g})$ and $\mathcal{S}_{0}$ was essential for understanding the negativity of the ${N_{\mathrm{c}}}$-terms. In the CXY approximation these terms turn into $\sim\alpha_{\mathrm{s}}C_{\mathrm{F}}(1/\varepsilon)\left(\mathcal{S}(x_{g})-\mathcal{S}_{0}\right)$, which can similarly be neglected at this level of accuracy. For the ${N_{\mathrm{c}}}$-terms the choice between $X(\xi)$ and $x_{g}$ was dictated by the desire to maintain the relation between the integral form of the BK equation and the subtraction. For the $C_{\mathrm{F}}$-terms there is no such requirement. Our choice here is to keep the same energy dependence $X(\xi)$, and the same kinematical limit $\xi<1-x_{g}/x_{0}$ for the $C_{\mathrm{F}}$-terms as for the ${N_{\mathrm{c}}}$-terms.

Because $\mathcal{I}$ and $\mathcal{I}_{v}$ vanish at $\xi=1$, replacing the upper limit $\xi<1-x_{g}/x_{0}$ by $\xi<1$ would be a rather good approximation here as long as we are safely inside the domain of validity of the eikonal approximation $x_{g}\ll x_{0}$. Evaluating the dipoles at different rapidity arguments $X(\xi)$ vs. $x_{g}$ makes a larger difference. However also here the difference becomes smaller at smaller $x_{g}$, because the integral over $\xi$ is increasingly dominated by smaller $\xi$ where $X(\xi)\approx x_{g}$. The extremely forward kinematical limit $x_{p}\to 1$ is an exception to this argument, because the vanishing of the collinear quark distribution at $x_{p}/\xi=1$ suppresses the small $\xi$ contribution to the cross section, making the result more sensitive to the difference between $X(\xi)$ and $x_{g}$. We will demonstrate the effect of these different choices for the $C_{\mathrm{F}}$-terms in a particular kinematical configuration in the next section.

III Numerical results: fixed coupling

We now turn to numerical results for the quark production cross section. In practice we consider up quark production and use the MSTW2008 NLO parametrization Martin:2009iq for the quark distribution $q(x)$ in the projectile. This distribution is evaluated at the transverse scale $Q=k_{\perp}$. The center of mass energy of the collision and the rapidity of the produced quark are chosen as $\sqrt{s}=500$ GeV and $y=3.2$ respectively. In these kinematics, the limit $x_{p}=1-x_{g}/x_{0}$ is reached when $k_{\perp}\sim 17.5$ GeV for $x_{0}=0.01$. As explained previously, we do not consider our calculation to be reliable for transverse momenta larger than this value so we will not show results above this limit.

We first consider a fixed value of the strong coupling, $\alpha_{\mathrm{s}}=0.2$. The rapidity dependence of the dipole operator of the target is obtained by solving numerically the LO BK equation with $\alpha_{\mathrm{s}}=0.2$ and an MV initial condition McLerran:1993ni ,

[TABLE]

with $Q_{\mathrm{s,0}}^{2}=0.2$ GeV2 and $\Lambda_{\mathrm{QCD}}=0.241$ GeV. Combining Eqs. (II.1) and (15), the NLO multiplicity reads

[TABLE]

As discussed above, recovering the LO BK result from this formulation is not a matter of dropping terms that are explicitly proportional to $\alpha_{\mathrm{s}}$. Instead this is done by setting $\xi=1$ everywhere in the integrand, except in the $\mathrm{d}\xi/(1-\xi)$ term that gives the leading logarithmic contribution, and in $X(\xi)$ which is needed to recover the BK equation in integral form. The separation between LO and NLO terms becomes more explicit if we use the equivalence between (II.1) and (13) to rewrite (19) as

[TABLE]

This last expression is very similar to the subtracted CXY cross section, except for the rapidity scale at which the dipole operators in the NLO corrections are evaluated and the upper limit on the $\xi$-integrals. We will see in the following that these changes lead to very different numerical results at large $k_{\perp}$ compared to the CXY expressions.

The fact that Eqs. (II.1) and (13) are equivalent is demonstrated explicitly in Fig. 1, where we show the LO+${N_{\mathrm{c}}}$ result obtained with these two equations. We observe that the results obtained with (II.1) are more stable numerically at large transverse momentum than when using (13). This was to be expected since (13) involves a cancellation between rather large contributions. One can also note that the NLO corrections proportional to ${N_{\mathrm{c}}}$ are rather large and negative at large $k_{\perp}$ but the LO+${N_{\mathrm{c}}}$ multiplicity is positive at all transverse momenta. This is in contrast with the CXY result, also shown on this figure, which becomes negative around $k_{\perp}\sim 6$ GeV.

Let us now consider the LO+$C_{\mathrm{F}}$ result,

[TABLE]

As already mentioned, the rapidity scale at which the dipole correlators are evaluated in the $C_{\mathrm{F}}$-terms is left unspecified by the proposal Iancu:2016vyg , and the difference between evaluating these correlators at $X(\xi)$ or $X(0)=x_{g}$ is formally a higher order effect. In Fig. 2 we show how the LO+$C_{\mathrm{F}}$ result is affected by these different choices. Note that at similar rapidities and transverse momenta at LHC energies $x_{g}$ (and $x_{p}$) would be smaller, and the difference between these choices smaller. In the following we choose to evaluate the dipole correlators of the $C_{\mathrm{F}}$-terms at $X(\xi)$ since it would be quite unnatural to use different rapidity scales in the $C_{\mathrm{F}}$ and ${N_{\mathrm{c}}}$-terms.

Now that we have specified how to evaluate both the ${N_{\mathrm{c}}}$ and $C_{\mathrm{F}}$-terms, we can sum these contributions to get the full NLO result. Since both the LO+${N_{\mathrm{c}}}$ result and the NLO $C_{\mathrm{F}}$ corrections are positive at all $k_{\perp}$, as shown in Figs. 1 and 2 respectively, the NLO multiplicity is positive as well. This is shown in Fig. 3, where one can observe that the net result of adding the ${N_{\mathrm{c}}}$ and $C_{\mathrm{F}}$ NLO corrections makes the cross section significantly smaller at high $k_{\perp}$ than the leading order result. For comparison we also show the corresponding CXY result which becomes negative above $k_{\perp}\sim 10$ GeV.

IV Numerical results: running coupling

A realistic QCD calculation must include the running of the coupling constant as a function of the scale. In the BK evolution equation this is typically done using the “Balitsky” prescription Balitsky:2006wa , as is the case for the dipole operators from Lappi:2013zma that we will be using here. Other choices are also possible (see e.g. Albacete:2004gw ; Kovchegov:2006vj ; Iancu:2015joa ), but the difference between these is not essential for our discussion here. What is important is that the BK equation is normally solved in coordinate space, and therefore also the running coupling prescription involves a kernel with a coupling constant that depends on coordinate differences. While momentum space solutions are also possible, there are several reasons why the coordinate space formulation is preferable. Firstly, the initial condition for the evolution needs to be extracted from a fit to experimental data, typically from DIS, where the cross section is most naturally expressed in terms of the coordinate space dipole. The fitting procedure is thus most naturally done in coordinate space. Secondly, a fundamental property of the BK equation is gluon saturation, which is guaranteed by the coordinate space dipole amplitude $S({\mathbf{r}})$ taking values between [math] and $1$. In a coordinate space calculation this unitarity requirement is explicitly satisfied also for running coupling. On the other hand, momentum space running coupling versions of the equation, at least known ones, are equivalent to a coordinate space version only parametrically, but not exactly. They do not therefore automatically enforce unitarity. Because of these reasons we will also here take the stand that the BK evolution that one would prefer to use (in our case that of Lappi:2013zma ) in conjunction with the cross section calculation will be a coordinate space running coupling one.

Parametrically the scale of the coupling constant is set in this case by the transverse momentum of the produced particle. The procedure suggested in Ref. Iancu:2016vyg is to simply replace the explicit coupling constant $\alpha_{\mathrm{s}}$ appearing in Eqs. (II.1), (13) and (15) by $\alpha_{\mathrm{s}}(k_{\perp}^{2})$. However, doing this while using dipole operators $\mathcal{S}$ whose energy dependence results from a BK equation with coordinate space running coupling spoils the exact equivalence between the manifestly positive (19) and subtracted (20) expressions for the cross section. This mismatch can be studied by comparing the LO limits of Eqs. (II.1) and (13). For this we use dipole correlators obtained by solving numerically the LO BK equation with the Balitsky prescription for the running coupling corrections Balitsky:2006wa . The initial condition is the MVe parametrization introduced in Lappi:2013zma :

[TABLE]

and the running coupling in coordinate space is taken as

[TABLE]

with $\beta_{0}=(11{N_{\mathrm{c}}}-2{n_{\mathrm{f}}})/3$ and $\Lambda_{\mathrm{QCD}}=0.241$ GeV. At large $|{\mathbf{r}}|$ the coupling is frozen at the value 0.7. The values of the parameters in these expressions were obtained in Lappi:2013zma by a fit to HERA DIS data Aaron:2009aa as $Q_{\mathrm{s,0}}^{2}=0.06$ GeV2, $C^{2}=7.2$ and $e_{c}=18.9$. This parametrization describes a target proton and not a nucleus, but this doesn’t affect our discussion. More importantly, the MVe parametrization has a positive Fourier transform.

In Fig. 4 we show the ratio between the LO limits of Eqs. (II.1) (involving a momentum space running coupling) and (13) (BK with “Balitsky”-coordinate space $\alpha_{\mathrm{s}}$). Here we use the expression

[TABLE]

for the momentum space running coupling that appears explicitly in the equations for the cross section, using different values of the parameter $C_{\text{mom}}^{2}$. From Fig. 4 we see that no value for $C_{\text{mom}}^{2}$ can give a good agreement between the two results for all $k_{\perp}$. In the following we fix $C_{\text{mom}}^{2}=10^{3}$ for which the disagreement is smaller than $30\%$ at both small and large $k_{\perp}$. This disagreement leads to significantly different results when evaluating the NLO cross section using either Eq. (19) or (20), as shown in Fig. 5. One immediately notices that the equivalence of the two equations is significantly broken, and the “subtracted” cross section (20) again becomes negative, albeit at a higher scale than in the CXY prescription. This demonstrates the point emphasized in Iancu:2016vyg that the subtraction procedure involves a difference of two terms that are numerically, although not parametrically, large and easily becomes unstable with any further approximations.

To solve this problem, one could in principle rewrite the expression for the cross section in coordinate space. This would allow to use a coordinate space running coupling matching the one used when solving the BK equation. However, how to implement the running coupling in a way which is consistent with the prescription Balitsky:2006wa is not unique, and a rather straightforward implementation leads to troublesome results (see the discussion in Appendix A). Therefore, we consider for now that the most physical prescription to use running coupling in this process is the unsubtracted result shown in Fig. 5.

V Outlook

In conclusion, we have demonstrated that it is possible to calculate single inclusive particle production in the hybrid formalism to NLO accuracy in a controlled way. While this is a promising start, many further improvements need to be made before a full phenomenological analysis and a comparison to experimental data can be made.

Firstly it is necessary to also include the other channels in the calculation. In principle this should be straightforward. In particular there is no reason why the same rapidity factorization procedure as here could not also be used for the $g\to g$ channel, which exhibits a similar large energy logarithm that needs to be resummed using the BK equation. One also needs to add fragmentation functions to the calculation. A successful phenomenology would require better control of also these, since so far there are large differences between different sets, in particular in LHC kinematics where the gluon channel starts to dominate d'Enterria:2013vba .

For a fully consistent NLO calculation one also needs to use a NLO version of the BK equation Balitsky:2008zza ; Lappi:2015fma ; Lappi:2016fmu or at least one including the resummed double and single transverse momentum logarithms Iancu:2015vea ; Iancu:2015joa that can be made to include most of the NLO effects Lappi:2016fmu . Using these requires the corresponding additional term to be added to the calculation of the single inclusive cross section. One also needs to obtain the initial condition for the BK equation from a fit to DIS data in an NLO calculation Balitsky:2010ze ; Beuf:2011xd ; Boussarie:2014lxa ; Beuf:2016wdz ; Boussarie:2016ogo . Also the issue of momentum vs. coordinate space running coupling merits to be studied further. A better consistency between calculations of DIS and particle production cross sections would require a consistent running coupling scheme between the two.

The NLO corrections almost inevitably decrease the cross section from the leading order calculation. This is essentially due to the fact that, as emphasized in Stasto:2014sea ; Ducloue:2016shw , at large $k_{\perp}$ the ${N_{\mathrm{c}}}$-part of the unsubtracted NLO corrections behaves like $\sim\xi/k_{\perp}^{4}$. The LO calculation only uses this gluon emission in the soft limit $\xi=1$ to construct the BK evolution equation, overestimating the result in exact kinematics with arbitrary $\xi$. Considering that the LO calculations such as Lappi:2013zma typically require “$K$-factors” greater than unity to make contact with the data, this might sound problematic. Here, however, it is important to keep in mind that also the NLO corrections to DIS cross sections will presumably be negative. Thus an NLO fit to DIS data could be expected to increase the normalization of the dipole amplitude compared to the current LO ones. An accurate assessment of the full phenomenological effect of the NLO corrections is impossible without a simultaneous calculation of DIS and forward pA particle production cross sections.

Acknowledgments

We thank E. Iancu and D. Zaslavsky for discussions and H. Mäntysaari for providing his BK solutions. This work has been supported by the Academy of Finland, projects 267321, 273464 and 303756 and by the European Research Council, grant ERC-2015-CoG-681707.

Appendix A Running coupling in coordinate space

The equivalence between the “unsubtracted” (19) and “subtracted” (20) expressions for the NLO multiplicity holds only if the coupling $\alpha_{\mathrm{s}}$ used to evaluate these expressions is the same as the one used when solving the Balitsky-Kovchegov equation to obtain the rapidity dependence of the dipole correlators. In particular, as shown in Sec. IV, this equivalence is broken if one wants to use a momentum space running coupling when evaluating the multiplicity while the prescription for the coupling used when solving the BK equation is formulated in coordinate space. In principle this problem could be solved by rewriting the expression for the multiplicity in coordinate space. A straightforward calculation shows that we can write the needed integrals $\mathcal{I}^{\text{finite}}$, $\mathcal{I}_{v}^{\text{finite}}$, $\mathcal{J}$ and $\mathcal{J}_{v}$ as

[TABLE]

With these coordinate space expressions it is possible to move the coupling inside the integrals and replace the kernels by an expression that (a) reduces to the same expressions for a fixed coupling and (b) reduces to a desired coordinate space running coupling BK kernel in the limit $\xi=1$. These two requirements do naturally not uniquely determine the expression at other values of $\xi$. In the notations introduced above the BK equation can be written as

[TABLE]

at the rapidity $Y=\ln\frac{1}{X}$. The Balitsky prescription Balitsky:2006wa for the running coupling consists in making the replacement $\alpha_{\mathrm{s}}\to\alpha_{\mathrm{s}}({\mathbf{r}}^{2})$ and in a modification of the kernel, leading to

[TABLE]

Therefore we find that a rather straightforward way to implement this prescription for arbitrary $\xi$ is to replace $\widetilde{\mathcal{J}}_{v}$ with

[TABLE]

as well as to move the coupling $\alpha_{\mathrm{s}}$ appearing in Eqs. (II.1), (13) and (15) inside the ${\mathbf{r}}$-integrals (A)-(A), then making the replacement $\alpha_{\mathrm{s}}\to\alpha_{\mathrm{s}}({\mathbf{r}}^{2})$. It is easy to check that the two conditions (a) and (b) mentioned previously are satisfied with this choice, which however is not unique. The numerical implementation of the “unsubtracted” (19) and “subtracted” (20) expressions for the NLO multiplicity using this coordinate space formulation is shown in Fig. 6. This figure demonstrates that, contrary to the results obtained with a momentum space running coupling shown in Fig. 5, the unsubtracted and subtracted results are the same, meaning that our modification of the ${N_{\mathrm{c}}}$-terms matches the Balitsky prescription at $\xi=1$. However, these results are problematic since the behavior of the NLO multiplicity is totally different from the results obtained at fixed coupling or with a momentum space running coupling, with a NLO result orders of magnitude larger than the LO one. This issue may be related to the fact that our implementation of the coordinate space running coupling, while quite straightforward, is not unique and other choices could give different results. How to consistently generalize the Balitsky prescription to $\xi\neq 1$ is a non-trivial issue that goes beyond the scope of this work.