Complete Next-to-Leading Order Calculation of Single Inclusive π 0 Production in Forward Proton-Nucleus Collisions

We present the first fully consistent calculation of inclusive π 0 production at forward rapidities at next-to-leading order (NLO) accuracy in proton-lead collisions at √ s = 8 . 16 TeV within the Color Glass Condensate approach. The center-of-mass energy dependence is determined by the Balitsky-Kovchegov equation with the initial condition constrained at NLO accuracy by the proton structure function measurements. We find that the LHCb data further constrain this non-perturbative input, and that both the NLO corrections to the evolution equation and to the impact factor have significant effects on the nuclear modification factor R p Pb . The complete NLO calculation is shown to be in qualitative agreement with the LHCb data, and the predictions are found to be sensitive to the details of the applied deep inelastic scattering (DIS) fit. These results demonstrate the need to perform global analyses of DIS and LHC data in order to probe gluon saturation phenomena at precision level at very small momentum fraction x .


I. INTRODUCTION
Proton-nucleus (pA) collisions at high center-of-mass energies studied at RHIC and especially at the LHC provide access to the part of nuclear wave function with very small momentum fraction x.This, together with the fact that densities in heavy nuclei are enhanced by a factor ∼ A 1/3 , makes forward particle production measurements in p+A collisions a promising channel to probe non-linear QCD dynamics and look for signatures of gluon saturation [1].
Inclusive particle production measurements at the LHC both at midrapidity [2,3] and at forward rapidities [4,5], together with the corresponding measurements at RHIC in forward kinematics [6][7][8], have shown clear modifications in the p+A cross section compared to the scaled proton-proton results.These modifications take place at transverse momenta of a few GeV, which is of the same order as the expected saturation scale for heavy nuclei at LHC energies.Thus, an accurate description of nuclear effects seen at the LHC is a crucial test for the gluon saturation picture at the highest collider energies available.To describe QCD in the high-density domain where saturation effects are important, an effective theory of QCD known as the Color Glass Condensate (CGC) [9,10] has been developed.It is the theoretical framework used in this paper, and leading order (LO) CGC calculations [11][12][13][14][15][16][17][18][19][20] have been successful in describing, at least qualitatively, different particle production spectra from RHIC and LHC.
Quantitative comparisons to the current and upcoming precision data require theoretical CGC calculations to be promoted to the NLO accuracy.In recent years, there has been a rapid progress in the field that has made NLO calculations in the saturation region possible.The crucial developments in the context of this work include the derivation of the NLO cross section for the single inclusive particle production obtained in Refs.[21,22], which are based on the leading-order formalism developed in [11,12].In the first phenomenological application, uncontrolled large NLO corrections were found rendering the cross section negative [23].This issue triggered a lot of attention in the community [24][25][26][27][28][29][30].To date, two solutions have been proposed.The first is through threshold resummation [28], which has led to an agreement with LHCb [4] and BRAHMS [7] data in Ref. [29] (however no NLO corrections to the high-energy evolution were included in [29]).The second solution to the negativity problem is by modifying the factorization procedure as shown in Ref. [25].This solution was extended to include the running coupling in [26].
In CGC calculations, the q q dipole-target scattering amplitude, a.k.a. the dipole amplitude, is a convenient degree of freedom used to describe interactions with the target at high energies.The center-of-mass energy dependence of the dipole amplitude, and that of the cross section, can be calculated perturbatively by solving the Balitsky-Kovchegov (BK) equation [31,32], which is also available at NLO accuracy [33][34][35].However, the starting point of the evolution equation requires a non-perturbative input describing the structure of the target proton or nucleus at moderate momentum fraction x.Typically parametrizations based on the McLerran-Venugopalan (MV) model [36][37][38] (which is strictly speaking valid for very large nuclei) have been used.The free parameters for the proton case have been extracted from fits to the HERA structure function data [39,40], and such fits at NLO accuracy became possible when the CGC results for the total deep inelastic scattering (DIS) cross section at NLO were obtained in Refs.[41][42][43].Recently, an alternative approach to obtain the initial condition for the BK evolution has been proposed in Refs.[44,45], where this initial condition is obtained based on the proton wave function in the valence region, including also the first perturbative correction from a gluon emission.
The purpose of this work is to perform the first fully consistent NLO calculation of inclusive π 0 production in LHC kinematics within the CGC framework.The initial condition for the small-x evolution at NLO is obtained from NLO fits to HERA data performed in Ref. [46].The determined dipole-proton amplitude is generalized to the dipole-nucleus case by employing the optical Glauber model [14], and as such there is no free parameter describing the nuclear high-energy structure despite the standard Woods-Saxon geometry.The π 0 production cross section is then calculated at the NLO accuracy by adopting the prescription for the factorization procedure and running coupling prescription developed in Refs.[25,26].
This work is structured as follows.In Sec.II, we present the calculation of inclusive π 0 production cross section at next-to-leading order within the CGC framework.The cross section is expressed in terms of the dipole-target scattering amplitude, whose initial condition and energy evolution are discussed in Sec.III.The results at the parton and hadron levels, together with comparisons to the LHCb data, are shown in Secs.IV and V, respectively.Finally, we present our conclusions in Sec.VI.A detailed discussion on the sensitivity of the results to the running coupling prescription in the hard impact factor is presented in Appendix.A.

II. PARTICLE PRODUCTION IN PROTON-NUCLEUS COLLISIONS
The forward particle production cross section in pA collisions at next-to-leading order is calculated in [21,22], where the cross section is shown to factorize to a hard part H a→c (a.k.a.impact factor), the part describing the eikonal interaction of partons with the target color field S a,c , and the collinearly factorized parton distribution function f a and fragmentation function D h,c : The process described in Eq. ( 1) involves a parton a from the incoming proton interacting with the nuclear target A, resulting in parton c that comes out and fragments into the observed hadron h.In this work, we include gluons and light quarks as we consider only π 0 production.The quark mass corrections are not included in the NLO impact factor [21,22] nor in the NLO DIS fit [46].Since the interaction with the nuclear target is eikonal, the shockwave picture applies.This allows S a,c to be expressed in terms of Wilson lines describing the eikonal propagation of partons in the target color field.In the large-N c limit employed throughout this work [47], each term of S a,c , involving Wilson line correlators, can be written in terms of the dipole-target scattering matrix which is a correlator of two fundamental light-cone Wilson lines V x ⊥ : where Here, P is a path ordering operator and A + is the classical color field, with the light-cone coordinates defined such that The proton is moving predominantly in the light-cone minus direction, while the target nucleus is moving in the plus direction.In Eq. ( 2), the angle brackets, ⟨• • • ⟩ X , denote the average taken over the target color field configurations.The correlator is evolved in Bjorken x from the initial X 0 down to X by solving the BK equation.
The hard impact factor, H a→c , includes perturbative corrections that can be calculated order-by-order.At the leading order, the incoming parton a and the fragmenting parton c are the same, a = c.This gives two possible channels of interaction: "quark channel" and "gluon channel", corresponding to whether a is a quark or a gluon.At NLO the quark can emit a gluon, or the gluon can split to a q q or gg pair.This leads to four possible channels at NLO: qq channel, qg channel, gq channel and gg channel, corresponding to a and c, respectively.The "primary parton" denoted by c that carries the momentum fraction ξ relative to a fragments into the hadron h.The phase space of the other unobserved parton is integrated over.The transverse integrals result in collinear divergences that are subtracted through the DGLAP evolution to the parton distribution functions (PDFs) and fragmentation functions (FFs).For consistency, we also use the PDFs and FFs at NLO.Specifically, throughout this work, we employ the MSTW PDFs [48] and the DSS fragmentation functions [49].
In the hard factor, there is also an integral over the longitudinal momentum fraction, ξ, of the primary parton.This integral contains rapidity divergent terms that can be included in the BK evolution of the target.In this work we apply the "unsubtracted scheme" as in Ref. [26] (and e.g. in related NLO calculations in Refs.[46,[50][51][52]), in which the rapidity divergence is not explicitly subtracted from the NLO contribution of H a→c .Specifically, dipole amplitudes in the LO term are evaluated at the initial scale, X 0 , of the small-x evolution.At NLO, dipoles are evaluated at where 0 < ξ < 1 − X g is the fraction of the primary parton longitudinal momentum carried by the fragmenting parton, and the fraction of the target longitudinal momentum transferred in the process is Here k ⊥ is the transverse momentum of the produced quark or gluon, and k ⊥ = |k ⊥ |.The kinematical upper bound ξ < 1 − X g ensures that the integral over ξ in the NLO impact factor does not diverge.This choice of limit is in accordance with the DIS fits [46] whose results are employed in the calculation of the dipole amplitudes in this work (see discussion in Sec.III).As our calculation is in the LHC kinematics where X g ∼ 10 −5 , the sensitivity to the large-x region is small.Before we proceed to the actual calculation, we introduce another important kinematic variable: the fraction of proton longitudinal momentum carried by the incoming parton before the emission of the primary parton.The ratio is given by In this scheme, at LO in the impact factor, the quarkand gluon-channel cross sections are given by [11,12,26,53] respectively.Here, the lower limit, τ = (p ⊥ / √ s)e y , is imposed to the fragmentation function integral so that x p ≤ 1, and p ⊥ = zk ⊥ .Note that Eqs.(7) are written in the large-N c limit.Furthermore, b ⊥ = (x ⊥ +y ⊥ )/2 is the impact parameter and r ⊥ = x ⊥ − y ⊥ .With proton targets (when calculating the cross section in proton-proton collisions), we assume, as in DIS fits [46], that the impact parameter dependence factorizes and replace Here, the proton transverse area σ 0 /2 is extracted in fits to proton structure function data.In the complete calculation that will sum over all channels, the light quark flavors in the quark-channel cross section (7a) will be summed over, accounting for the fact that a quark (or an antiquark) of any flavor inside the proton can interact with the heavy nuclear target.
As for the NLO in the impact factor, the cross section for the qq channel in the unsubtracted scheme has been calculated in Ref. [26] based on Ref. [21].The explicit expression for the cross section reads where We perform a similar calculation to obtain the cross sec-tion for the other three remaining channels (qg, gq, gg) in the unsubtracted scheme starting from the results ob-tained in Ref. [21].The resulting expressions are where Similarly to the LO case, in the actual cross-section calculation, every channel except for the gg channel involves a sum over the light quark flavors.
In Eqs. ( 9) and ( 11), the strong coupling α s is taken to run with the transverse momentum of the interacting parton throughout this work.This scheme choice was shown in Ref. [26] to avoid the "fake potential" issue that would result in uncontrolled NLO correction.In Appendix A, we demonstrate that there is only a moderate sensitivity on the running coupling scheme that mostly cancels in cross section ratios, as long as we use one of the schemes that avoid the fake potential issue.
The z-integrals in Eqs. ( 7), (9), and ( 11) can be interpreted as the convolutions between the fragmenta-tion function and the parton-level cross section, which is the single inclusive parton production cross section in pA collisions.In some cases of the calculation below, the parton-level cross sections will be negative at given parton transverse momentum, k ⊥ .This negativity problem results from the fact that already at leading order the cross section is proportional to the Fourier transform of the dipole amplitude (see Eq. ( 7)), and this Fourier transform is not guaranteed to be positive definite. 1We will return to this discussion in Sec.III.To avoid propagating the negative cross section problem to the hadron level, we replace the parton-level cross section, dσ parton , by max{0, dσ parton } in the integrand of all the fragmentation-function integrals over z that appear in Eqs. ( 7), (9), and (11).We have confirmed numerically that this has a negligible effect on the hadron-level cross section except at p ⊥ values where the hadron-level results would be negative.

III. DIPOLE AMPLITUDE AND THE HIGH-ENERGY EVOLUTION
As we see from Eqs. ( 9) to ( 12), each term in the NLO cross section involves dipole amplitudes, S(r ⊥ , X(ξ)), for some transverse separation r ⊥ .In the shockwave picture, these amplitudes correspond to the interactions between the color-singlet quark-antiquark dipole and the target nucleus.The dipoles are evolved via the small-x BK evolution equation from the initial value, X 0 , to the smaller value of X(ξ) given in Eq. ( 4).
In the case where the nuclear target is in fact a proton, we adopt the parametrization employed in [46] (and previously in leading order fits [14,54]), where we take the non-perturbative initial conditions for the dipole-proton scattering amplitude to be Here Λ QCD is taken to be 0.241 GeV, while γ and Q 2 s,0 are free parameters that are determined by performing a fit to the proton structure function data from HERA [39,40].The determined dipole-proton scattering amplitudes are available at [46,55].Physically, Q 2 s,0 controls the saturation scale at the initial condition and γ is the anomalous dimension which determines the shape of the dipole amplitude in the limit of small dipole size r ⊥ .In addition, other free parameters include the proton transverse area, σ 0 /2, and a parameter controlling the scale of the running coupling in coordinate space.The latter is denoted by C 2 in Refs.[14,46,54], see Eq. (A1) for an explicit expression for the α s as a function of transverse distance scale.
The first fit to determine the required non-perturbative input at next-to-leading order accuracy was performed in Ref. [46], and the determined dipole amplitudes are also used in this work.In principle, in the fully consistent NLO analysis, one should also use the NLO BK evolution equation [33].As the NLO BK equation is computationally challenging, it was approximated in Ref. [46] by the leading order BK equation together with the running coupling corrections and a resummation of the most important higher order corrections enhanced by large transverse logarithms.This imposes a dependence on the resummation scheme that specifies how the large transverse logarithms that appear in the collinear regime of the NLO corrections are resummed.It is worth noting that this collinear resummation scheme employed to approximate the NLO BK equation is distinct and separate from the "subtraction scheme" of the cross section, which is discussed in Sec.II and Ref. [26].The three different resummation schemes, which lead to three different approximations of the NLO BK evolution considered here, are: the kinematically constraint BK (KCBK) equation [56], the rapidity local resummed BK (ResumBK) equation [57,58], and the target-momentum BK (TBK) equation [59].The TBK equation is formulated in term of the target rapidity, η, which is then transformed to the projectile rapidity used to derive the impact factor as in Ref. [46].Here, the terminology is the same as in Ref. [46], and we note that in Ref. [60] it was found that the remaining NLO corrections not enhanced by these large transverse logarithms have very small effects on DIS cross sections.These evolution equations have been successfully employed in NLO calculations of exclusive vector meson [50,51] and heavy quark [52] production in DIS.
In Ref. [46], many different fits corresponding to different starting points of the BK evolution and different running coupling schemes are used.Although the running coupling prescriptions employed in [46] do not exactly correspond to the ones used in the impact factor, the effect of the discrepancy is formally of higher orders in α s .Additionally, the reduced cross section data were fitted separately from the pseudodata consisting of an estimated light quark contribution.In this work, we only consider fits to the reduced cross section data, as all fits to the light quark pseudodata have a large anomalous dimension, γ > 1.This could make the cross section negative, as discussed later in this Section.Of the fits to the reduced cross section data, we use the ones for which the evolution starts at X 0 = 0.01 (above which the dipole is frozen at its initial condition), unless otherwise specified.
For the dipole-nucleus scattering process, it is not possible to perform a similar fit to the structure function data to obtain the scattering amplitude due to the lack of small-x nuclear-DIS data.Instead, we generalize the dipole-proton amplitude to the dipole-nucleus case using the optical Glauber model as in Ref. [14].This corresponds to modifying the initial saturation scale from Eq. ( 13) to account for the impact-parameter profile of the heavy nucleus.The dipole-nucleus amplitude that we use as an initial condition for the BK evolution then reads The free parameters (σ 0 /2, Q 2 s,0 and γ) are exactly the same as in the dipole-proton amplitude, and thus there is no freedom when moving to nuclear targets.Note that in the dilute regime the dipole-nucleus cross section is exactly A times the dipole-proton one, and S → 0 in the saturation regime.Here, A is the mass number of the target, and T A (b ⊥ ) is the transverse thickness function of the nucleus, which can be obtained from the Woods-Saxon distribution of nuclear density.It is normalized such that d 2 b ⊥ T A (b ⊥ ) = 1.The small-x evolution is then applied to the initial condition ( 14) separately for each b ⊥ , resulting in different dipole-nucleus interaction amplitudes as functions of x and r ⊥ for different values of b ⊥ .This method generalizing the dipole-proton amplitudes into their dipole-nucleus counterparts has been employed in [61] for vector meson production at NLO in γA collisions.Furthermore, a successful description of forward particle production cross section in protonnucleus collisions have been obtained in leading order calculations in Refs.[14,[16][17][18].
Another important characteristic of the fit in [46] is the fact that the anomalous dimension, γ, is allowed to exceed one.This is a similar practice to that employed in the preceding LO fits [14,54].Unlike the total DIS cross section that only depends on the dipole amplitude convoluted with the photon wave function, the inclusive parton production cross section is proportional to the two-dimensional Fourier transform of the dipole amplitude.This Fourier transform is highly sensitive to the detailed shape of the dipole amplitude at very small transverse dipole size, r ⊥ , which is the region that contributes minimally to the total DIS cross section.In particular, with γ > 1, the Fourier transform is not positive definite [62], and as a result the parton level cross section can be negative.In order to minimize this issue, we mostly focus on fits reported in Ref. [46] for which the anomalous dimension is γ ≈ 1.An additional advantage of fits that have γ ≈ 1 is the fact that only in this case does the nuclear saturation scale obey the natural scaling, Q 2 s,0,A = σ0 2 AT A (b ⊥ )Q 2 s,0 .Furthermore, as mentioned at the end of Section II, in the regions of parton transverse momentum, k ⊥ , where the parton-level cross section is still negative, we replace these negative cross sections by zero when evaluating the convolution with fragmentation functions.

IV. NUCLEAR MODIFICATION AT PARTON LEVEL
In order to determine how the NLO corrections in the evolution equation and in the impact factor affect the single inclusive cross section, we first calculate inclusive parton production cross section in the LHC kinematics.This means that the fragmentation function integral is not included.At this level, the gg channel dominates in the spectra, followed by the qq channel.While we will include all channels when calculating the hadron-level results in Sec.V, the parton level results are shown only for these two channels to illustrate the importance of the various NLO corrections.In this Section, we also focus on the most central collisions only, i.e. b ⊥ = 0 in the applied optical Glauber model, as in this case the nuclear saturation scale is the highest and nuclear effects are most pronounced.However, we note that the b ⊥ = 0 collisions do not correspond to the experimental definition of most central pA collisions.In Sec.V, we will present results at hadron level for minimum bias collisions in which case the applied optical Glauber model (integrated over all b ⊥ ) is expected to be applicable.As the genuine parton-level cross section is not finite, we calculate (as in Ref. [26]) only the part that does not involve the 1/ε divergence cancelled by the renormalization of the fragmentation function.
This exercise is motivated by the fact that the (approximative) next-to-leading order BK evolution employed in this work results in a different shape of the dipoletarget scattering amplitude compared to the leading order evolution.More precisely, the leading order BK evolution drives the solution towards an asymptotic shape N = 1 − S ∼ r 2γ ⊥ with γ ∼ 0.6 . . .0.8 at small dipoles.This evolution removes the Cronin enhancement at moderate transverse momenta that appears in leading order calculations when an MV model initial condition is used [63].On the other hand, the NLO evolution does not significantly modify the anomalous dimension, γ, i.e. the shape of the dipole amplitude in the very small r ⊥ region remains relatively close to that of the initial condition [34,[57][58][59].
The nuclear modification factors shown in Fig. 1 for the quark and gluon production at most central collisions (b ⊥ = 0) are defined as Here, is the number of binary nucleon-nucleon collisions computed from the optical Glauber model with σ inel referring to the total inelastic nucleon-nucleon cross section.The invariant yield in proton-nucleus collisions, dN pA qq,gg , is obtained by not performing the d 2 b ⊥ integral in expressions (9) and (11) shown in Sec.II and evaluating the dipole amplitudes at fixed b ⊥ = 0. On the other hand, the proton-proton yield is obtained from the proton-proton cross section as [14] dN pp qq,gg = 1 σ inel dσ pp qq,gg .
With all these ingredients, the quantity in Eq. ( 15) compares the pA yield to its pp counterpart and becomes 1 when the pA collisions can be seen as an incoherent superposition of independent proton-proton collisions.Note that the dependence on σ inel cancels in Eq. (15).To demonstrate the effect of NLO impact factor, we compare the full NLO calculation to the case where the LO impact factor is used with the NLO dipole.The latter case corresponds to Eq. ( 7) with the dipole amplitude evolved to X g .In both cases, the factorization scale for the PDFs is set to µ = 4k ⊥ , where k ⊥ is again the transverse momentum of the produced parton.We use the KCBK resummation scheme for the BK evolution with the parent dipole running coupling prescription.Since the NLO evolution does not evolve the dipole towards an asymptotic shape with a small anomalous dimension γ, we obtain a large Cronin enhancement peaked at k ⊥ ≈ 12 GeV in the forward LHC kinematics when the hard impact factor is calculated at the leading order.As discussed above, with a leading-order evolution one does not obtain a Cronin peak in the LHC kinematics [14,63].Hence, we conclude that NLO corrections to the small-x evolution (not included in previous NLO calculations for forward particle production [26,29]) can have a qualitatively significant effect on the nuclear modification factor.This Cronin peak vanishes when the NLO impact factor is used together with the NLO evolved dipole, highlighting the importance of performing the calculation consistently at the NLO accuracy.We note that no Cronin peak is seen in the forward LHCb [4,5] or midrapidity ALICE [2,3] data.
Furthermore, we find a strong nuclear suppression at low parton transverse momenta.The NLO results have a significantly weaker dependence on the gluon k ⊥ compared to the case where the leading order impact factor is used.This is because the latter corresponds to the 1 → 1 kinematics, in contrast to the 1 → 2 kinematics resulting from a hard parton emission vertex at NLO.In the NLO case, the unobserved parton carries a significant transverse momentum, rendering the k ⊥ dependence weaker, see also related discussion in Ref. [64].Finally, in the dilute (high-k ⊥ ) region, the nuclear effects vanish and R pPb → 1.This is in contrast to the leading order calculations where the leading order BK evolution changes the asymptotic shape of the dipole amplitude, and consequently one obtains R pPb → 1 at moderate transverse momenta k ⊥ ∼ 5 . . . 10 GeV only at the initial condition [14,63].At the LHC kinematics, a significantly large transverse momentum is needed in order to be sensitive to dipoles close to the initial condition where X g ≈ X 0 , see Eq. ( 5), and as such R pPb would approach unity much more slowly in a complete LO calculation.TeV and y = 3, as a function of the gluon transverse momentum.The KCBK resummation scheme is employed to approximate the complete NLO BK evolution, together with the parent-dipole running coupling prescription.The labels, LO and NLO, refer to the order in the impact factor.

V. HADRON LEVEL RESULTS
In this section, we calculate the complete forward π 0 production cross sections for minimum bias proton-lead collisions at √ s = 8.16 TeV.When calculating the impact parameter integral (recall that at each b ⊥ we use independently evolved dipole-nucleus amplitudes), we encounter regions where the saturation scale of the nucleus would be smaller than that of the proton.To avoid unphysically rapid growth of the nuclear size in this region we calculate the yield in pA collisions at large b ⊥ following Ref.[14] as This corresponds to R pPb = 1 at large impact parameters.
The inclusive π 0 production cross sections at NLO as functions of pion transverse momentum compared to the LHCb data [5] are shown in Figs. 2 separately for the different resummation schemes used to approximate the full NLO BK evolution (KCBK, ResumBK and TBK).To demonstrate the sensitivity to the running coupling prescription employed in the BK evolution and to the related effects of the parameters in the initial condition, we also show the KCBK results using two different running coupling schemes employed in the NLO DIS fits in Ref. [46].Note that in the NLO impact factor, Eqs. ( 9) and (11), the strong coupling is always evaluated at the scale given by the transverse momentum of the produced parton.The dependence on the factorization scale is illustrated by varying µ in the range µ = 2p ⊥ , . . ., 8p ⊥ .
The LHCb data shown for comparison span the rapidity range of 2.5 < y < 3.5.We calculate the cross section at the central value, y = 3. Generically, the p ⊥ dependence of data is well captured in all the employed setups, except when the "Balitsky+smallest dipole" running cou- pling scheme (as defined in Ref. [46]) is used in which case the spectrum falls more steeply than the data.Similarly, the overall normalization of the LHCb data is typically slightly overestimated, except in the Balitsky+smallest dipole running coupling setup in which case the normalization is also underestimated by a factor ∼ 2 . . .3. As we show in Appendix.A, the overall normalization depends slightly (up to ∼ 50%) on the running coupling prescription applied in the impact factor, but the shape of the p ⊥ spectra and the nuclear modification factor are insensitive to this scheme choice.We recall that in Ref. [46] all resummation schemes for the BK equation applied here result in identical γ * p cross sections.This demonstrates that single inclusive hadron production at the LHC provides additional complementary constraints to the extraction of the non-perturbative initial condition for the dipole-proton scattering amplitude.Similar conclusions have been made previously in studies involving exclusive vector meson production [50,51] and heavy quark production [52] in DIS.This complementarity is due to the fact that the π 0 cross section in forward pA collisions is sensitive to the dipole amplitude at different length scales than the total DIS cross section.In particular, small dipoles do not contribute to the structure function F 2 , but on the other hand the shape of the dipole amplitude at small dipole size can have a dramatic effect, even rendering the (parton-level) cross section negative if the dipole amplitude vanishes faster than r 2 ⊥ at small r ⊥ as discussed in Sec.III and in Ref. [62].
We note that in Ref. [46] more fits than the four that we use in Figs. 2 are reported.We will further demonstrate the sensitivity of our results to the details of the NLO DIS fit later in this Section, but note that only the fits with parent dipole prescription used in Figs. 2 have anomalous dimensions, γ ≈ 1.This is preferred as the parton level cross section is then positive definite and the optical Glauber model extension to the dipole-nucleus scattering is more natural.For comparison, we show in Fig. 2b the result obtained using the KCBK evolution with the Balitsky+smallest dipole running coupling prescription, a fit that has the anomalous dimension of γ = 1.21 at the initial condition.As discussed above, in the NLO evolution the anomalous dimension γ does not change significantly [34,46], and as such there remains a sensitivity to the initial anomalous dimension even at the LHC energies.Larger anomalous dimensions typically result in more steeply falling spectra as observed e.g. in Ref. [14] in a leading-order calculation.This is also evident in our NLO results when comparing Figs.2a  and 2b: the spectra in Fig. 2b with γ = 1.21 at the initial condition fall more steeply with p ⊥ than the LHCb data, compared to the case shown in Fig. 2a whose spectra describe the p ⊥ dependence of the LHCb measurement well with the initial anomalous dimension of γ ≈ 1.Other fits with even larger γ would yield p ⊥ spectra that fall even more steeply than the one shown in Fig. 2b, and hence these fits are not favored by the LHCb data.We note that, as shown in Ref. [14], the leading order calculations require a larger γ ∼ 1.2 in order to reproduce the p ⊥ dependence seen in the midrapidity proton-(anti)proton measurements at the LHC [65,66] and at Tevatron [67].
Let us next consider the nuclear modification factor at the hadron level, defined as This was studied in the case of most central collisions (b ⊥ = 0) at the parton level in Sec.IV.The complete hadron-level results are shown in Figs. 3 together with the corresponding LHCb measurements [5].The obtained nuclear modification factor is similar with the four NLO DIS fits applied.We find weak nuclear suppression at low hadron p ⊥ and obtain R pPb → 1 at moderate p ⊥ , unlike in leading order calculations [14].The amount of nuclear suppression is significantly underestimated in all applied setups.The origin of this can be traced back to the fact that the applied NLO fits for the dipole amplitude with γ ∼ 1 (required to get positive definite parton level cross section) prefer relatively small proton transverse area σ 0 /2 = 10 mb.This is significantly smaller than σ 0 /2 = 17 mb obtained from the leading-order fits [14,54].Since the saturation scale of the nucleus scales as Q 2 s,0,A ∼ σ0 2 AT A (b ⊥ )Q 2 s,0 according to the applied optical Glauber model (with γ = 1), c.f. Eq. ( 14), we have a smaller difference between the proton and nuclear saturation scales compared to the leadingorder case.This results in a weaker nuclear suppression at NLO.Furthermore, note that the applied fits with smaller σ 0 can be seen to be preferred by measurements of several other observables.For example, the momentum transfer dependence of the exclusive J/ψ production measurement from HERA [68,69] corresponds to σ 0 /2 = 10 mb [70], although the exact value depends on the assumed proton shape that cannot be exactly determined from the currently available data.Along the same line, the diffractive structure function data have been shown in Ref. [71] to require a very steep proton density profile when σ 0 /2 = 17 mb is used in leading-order calculations, suggesting that somewhat smaller proton size could be preferred.
In order to further demonstrate the sensitivity to the initial condition used in the small-x evolution, we calculate the inclusive π 0 production cross section and the nuclear suppression factor using an initial condition parametrization that has a larger proton transverse area, σ 0 /2 = 18.39 mb.This is the KCBK initial condition with parent dipole running coupling prescription where the BK evolution starts at X 0 = 1 in Ref. [46].The resulting pion spectrum and nuclear modification factor are shown in Figs. 4. Despite the relatively small difference one might expect from varying the starting point of the high-energy evolution, Figs. 4 display clear differences when compared to the respective counterparts in Figs. 2  and 3.In large part, this is due to the differences in the BK initial condition's parameters that resulted from the fit to the HERA data in Ref. [46].In the large-σ 0 case of Figs. 4, we find that the overall cross section is typically overestimated by a factor ∼ 2 as seen in the transverse momentum spectra shown in Fig. 4a.This initial condition has an anomalous dimension γ = 1.21, but unlike above when applied to fits where the BK evolution starts at X 0 = 0.01, the p ⊥ slope is only slightly too steep compared to the LHCb data.The fact that a larger γ is not clearly disfavored in this case can be traced back to the fact that there is now an additional BK evolution from X = 1 to X = 0.01, in contrast to the previous case where the BK evolution only starts at x = 0.01.Although the shape of the dipole amplitude does not change at asymptotically small r ⊥ , the evolution can have a moderate effect on the shape of the dipole at intermediate dipole sizes [46], which can result in slight modifications on the p ⊥ slope.
The nuclear suppression factor shown in Fig. 4b shows a much stronger suppression than those in Figs. 3, where fits with significantly smaller proton transverse areas, σ 0 /2 = 10 mb, were used.As already discussed above, with larger σ 0 /2 the relative difference between the proton and nuclear saturation scales is larger, resulting in stronger nuclear suppression.With this particular fit having σ 0 /2 = 18.39 mb, we obtain an even stronger nuclear suppression at low pion transverse momentum than the level observed in the LHCb data.This demonstrates that R pPb measurements are very sensitive to the proton size parameter σ 0 /2 (at least in the applied Optical Glauber model for the nucleus).This is in contrast to the DIS structure function case that are fitted to extract the initial condition for the BK evolution: in the linear regime the structure functions depend only on the product σ0 2 Q 2 s,0 and as a result significant correlations between these parameters can be expected.However, as no uncertainty estimate for the NLO DIS fits is currently available, it is not possible to quantify how the uncertainties in the fit parameters actually propagate to the nuclear modification factor calculated here.
Finally, we study how the cross section and the nuclear modification factor depend on the pion rapidity which controls the amount of BK evolution in the calculation.These results can be seen as predictions for the future measurements.The LHCb collaboration has previously measured inclusive charged hadron production up to y = 4.3 [4].Furthermore, in the future, the ALICE collaboration with the Forward Calorimeter (FoCal) detector upgrade should be able to perform measurements almost up to y = 6 [72].FIG. 4. Inclusive π 0 spectra and nuclear modification factor computed using the KCBK resummation scheme for the BK evolution and a fit where the BK evolution starts at X0 = 1.This particular fit has a larger proton transverse area, σ0/2 = 18.39 mb.The results are compared with the LHCb data [5].The scale variation band is similarly generated with µ ∈ {2p ⊥ , 4p ⊥ , 8p ⊥ }.
The inclusive π 0 spectra at rapidities y = 3, 4, 5, 6 are shown in Fig. 5, and similarly the nuclear modification factors at the same values of rapidity are shown in Fig. 6.Here, we only use the KCBK resummation scheme for the BK evolution with the parent dipole running coupling prescription in the case where the evolution starts at X 0 = 0.01, i.e. with the proton size σ 0 /2 = 9.74 mb and γ = 0.98, that was used also in Figs.2a and 3a.A reason for this particular choice of resummation scheme, running coupling prescription and evolution starting point lies in the fact that this combination results in γ < 1 according to the fit in Ref. [46].This implies that the cross section is positive definite as long as the unsubtracted scheme is employed in conjunction with the momentum-space run- ning coupling prescription for the hard factor [26].As discussed above, the magnitude of the nuclear suppression obtained depends strongly on the proton transverse area on which there is significant variation between the available DIS fits, but the rapidity dependence should not depend strongly on the proton size.The π 0 production cross section is heavily suppressed at large rapidities as the quark and gluon PDFs vanish rapidly when approaching the kinematical boundary x p = 1.The nuclear suppression is found to be stronger for more forward rapidities due to the larger saturation scale of the nucleus probed in such kinematics.In the moderately large transverse momentum region, p ⊥ ≳ 5 GeV, we again get R pPb → 1.This feature is independent of the pion rapidity and qualitatively different from what is seen in leading-order calculations [14].These observations qualitatively agree with the trends found in the LHCb measurements of charged hadron production at different rapidities (with √ s = 5 TeV) [4], although the rapidity dependence obtained in this work is weaker.

VI. CONCLUSION
In this work, we have performed the first calculation of single inclusive π 0 production completely at nextto-leading accuracy in the Color Glass Condensate approach consistently with the DIS structure function data.Forward particle production processes in p+A collisions probe nuclear wave functions in the small-x regime.As a result, significant saturation effects can be expected in this kinematics.The developments presented in this work enable precision studies of gluon saturation at the smallest values of x accessible with the currently available collider energies.
First, we extended the results of Ref. [26] and calculated contributions from all partonic channels to the inclusive particle production cross section in the unsubtracted scheme.This allowed us to consistently use the dipole-target scattering amplitudes extracted from the HERA structure function data at NLO in Ref. [46].These results are shown in Eqs.(11).Then, we demonstrated that the NLO corrections both to the small-x evolution and to the impact factor have qualitatively significant effects on the nuclear modification factor in inclusive particle production processes.A distinct qualitative feature of the fully consistent NLO results presented in this work is that, at all rapidities, the R pPb approaches 1 already at moderate hadron transverse momentum, p ⊥ = 5 . . . 10 GeV.This is in contrast to the LO results for which the R pPb remains suppressed up to very large p ⊥ in the LHC kinematics [14].
When both the impact factor and the small-x evolution are consistently taken at NLO accuracy, we obtain a nuclear modification factor and inclusive π 0 spectra that qualitatively agree with the LHCb data as shown in Figs. 2, 3 and 4. Depending on the details of the applied DIS fit parametrizing the initial condition for the highenergy evolution, in particular on the extracted proton (gluonic) transverse area, our results for the nuclear modification factor R pPb either under-or overestimate the observed nuclear suppression in the inclusive forward π 0 production.Overall, the shape of the transverse momentum spectra is found to further constrain the DIS fits, disfavoring parametrizations with large anomalous dimensions γ corresponding to steeply falling dipole-proton amplitudes at small dipole sizes.We note that similar additional constraints to the non-perturbative initial condition have recently been obtained when other DIS data sets have been included in the analysis in addition to the total cross section [50][51][52].
Because of the strong sensitivity to the initial condition of the BK evolution it is necessary to perform global analyses in the CGC framework that include a variety of different DIS and p+A observables with rigorously propagated uncertainties, in order to obtain robust signals of non-linear QCD dynamics at currently available collider energies.The feasibility of such analyses simultaneously including both the total DIS cross section and forward LHC data has been demonstrated in this work where the first consistent NLO calculation has been presented, and the importance of a consistent treatment of different NLO corrections has been illustrated.
In this Appendix, we compare the results obtained using different running coupling prescriptions in the NLO impact factor, Eqs. ( 9) and (11), while using the same parent-dipole prescription in the small-x evolution.For simplicity we only consider π 0 production in the most central b ⊥ = 0 case.Specifically, the running coupling prescriptions considered here are: (i) Fixed coupling, with ᾱs = α s N c /π = 0.2.
(ii) Parent-dipole prescription, given by the transverse separation between the positions of the incoming partons (before the primary parton emission) in the amplitude and the complex-conjugate amplitude.The coupling constant in term of a transverse separation, r ⊥ , is given by [46] α s (r ⊥ ) = 4π where β 0 = (11N c − 2N f )/3 and Λ QCD = 0.241 GeV.Here, recall that we include only the three lightest quark flavors: N f = 3.For the parameters characterizing the infrared behavior of the coupling constant, we follow the convention from [46] and choose µ 0 /Λ QCD = 2.5 and c = 0.2.
(iii) Momentum-space prescription, given by the transverse momentum, k ⊥ , of the produced parton, such that α s (k ⊥ ) = 4π with the same values for β 0 and Λ QCD .
(iv) Mixed-space prescription, which is identical to the momentum-space prescription, except for the qq-channel terms proportional to J and J v , c.f. Eqs. ( 9), (10a) and (10b).In those terms, the coupling constant runs with the separation, x ⊥ , between the two partons after the primary parton emission, according to Eq. (A1).This is the prescription employed in [26].In this work, we find that a finite result can only be achieved when we restrict position-space prescriptions to these two specific J and J v terms; in other terms, any linear combination of transverse positions involving the daughter dipole would contain an undesirable "fake potential" contribution, which will be discussed shortly.Hence, in most terms, we resort to the momentum-space prescription given in Eq. (A2).
The π 0 spectra for pPb collisions are plotted in Fig. 7 using the four different running coupling prescriptions in A plot of pPb spectra computed using different running coupling prescriptions in the impact factor, while using the parent-dipole prescription in the BK evolution.
the impact factor.There, we see that the momentumspace and mixed-space prescriptions result in spectra that differ from each other by an overall normalization factor, but the shape is almost identical.On the other hand, when the running coupling effects are neglected in the fixed coupling case, one obtains a somewhat harder spectrum as expected.Overall differences in the cross sections calculated using these three schemes are within 50%.In contrast, a stark difference is observed when the parent-dipole prescription is used.In that case, the spectrum is much greater in magnitude and displays an exponential decay with p ⊥ (note the logarithmic scale in Fig. 7).The effect is also clearly visible in the nuclear modification factor.Fig. 8 shows the R pPb for different running coupling prescriptions in the impact factor.Here, we see that the results are almost identical for all the prescriptions except for the parent-dipole, whose nuclear modification factor simply remains a constant.
The qualitatively different behavior obtained with the parent-dipole scheme is due to the "fake potential" prob-

Momentum-space Fixed
Mixed-space Parent-dipole FIG. 8.A plot of nuclear modification factor, R pPb , computed using different running coupling prescriptions in the impact factor, while using the parent-dipole prescription in the BK evolution.Note that the blue and magenta lines, corresponding respectively to the momentum-space and fixed coupling prescriptions, lie on top of each other.lem discovered in Ref. [26].The main idea is that the dependence of the coupling constant on a transverse position changes the Fourier transform of the dipole amplitude, creating extra unphysical terms that dominate the cross section.The constant in the parent-dipole R pPb is simply the ratio between the unphysical terms that would result from the pPb and pp initial conditions.Finally, it is worth noting that employing any positionspace prescription in any of the terms besides the J and J v terms in the qq channel discussed above would also produce a fake-potential contribution, resulting in the spectra and nuclear modification factor identical to those of the parent-dipole prescription.
With the relatively minor differences among all the prescriptions that are free of the fake-potential contribution, we decide to employ the momentum-space prescription for the impact factor throughout the main sections of this work.

FIG. 1 .
FIG.1.Nuclear modification factor in the dominant qq and gg channels for the most central proton-lead collisions calculated at the LHCb kinematics, √ s = 8.16 TeV and y = 3, as a function of the gluon transverse momentum.The KCBK resummation scheme is employed to approximate the complete NLO BK evolution, together with the parent-dipole running coupling prescription.The labels, LO and NLO, refer to the order in the impact factor.

FIG. 2 .
FIG.2.Inclusive p + Pb → π 0 + X cross section as a function of pion transverse momentum, p ⊥ , using different resummation schemes and different running coupling prescriptions in the BK evolution compared with the LHCb data from Ref.[5].Here, √ s = 8.16 TeV and the pion rapidity is y = 3.The blue band is obtained by varying the factorization scale from µ = 2p ⊥ to µ = 8p ⊥ , and the central lines is obtained with µ = 4p ⊥ .The statistical and systematic uncertainties for the LHCb data are added in quadrature.

FIG. 3 .
FIG.3.The nuclear modification factor R pPb obtained using different resummation schemes and running coupling prescriptions in the BK evolution compared to the LHCb data[5].The results are calculated at √ s = 8.16 TeV and y = 3.In each plot, the blue band with the black center line represents our numerical computation with µ ∈ {2p ⊥ , 4p ⊥ , 8p ⊥ }.The statistical and systematic uncertainties for the LHCb data are added in quadrature.
FIG. 7.A plot of pPb spectra computed using different running coupling prescriptions in the impact factor, while using the parent-dipole prescription in the BK evolution.