Probing parton saturation with forward Z-boson production at small transverse momentum in p+p and p+A collisions

We calculate and compare the differential cross sections for forward Z-boson production at small transverse momentum, in proton-proton and proton-nucleus collisions, using both the collinear and dilute-dense factorization frameworks. In both cases, we implement a Sudakov resummation of the large logarithms generated by soft-gluon emissions, which is essential in order to describe the transverse momentum distribution of forward Z bosons measured at the Tevatron and the LHC. We further compute the nuclear modification factor in the dilute-dense framework, hoping to single out signals of saturation effects at small values of x. Our predictions are compared with those obtained in the collinear factorization framework, using two different nuclear parton distribution functions.


I. INTRODUCTION
The Large Hadron Collider (LHC) offers unique opportunities to study the QCD dynamics at low partonic longitudinal momentum fractions x. Of particular interest are observables involving particle production in the forward region, which allow to probe parton densities in one of the colliding hadrons down to x ∼ 10 −5 . On theoretical grounds, at such low values, one expects small-x effects to be relevant, in particular the high-gluon-density phenomenon known as the gluon saturation [1][2][3][4]. These asymmetric processes can be described in a hybrid approach [5][6][7], in which the well-understood projectile, made of a dilute parton content, is described in terms of standard collinear parton distribution functions (PDFs), while the small-x dynamics in the saturated target wave function is dealt with using the Color Glass Condensate (CGC) effective theory. Then, comparing proton-proton (p+p) collisions and proton-nucleus (p+A) collisions, one can predict systematic features of non-linear QCD dynamics at high parton densities, and look for them in experiments.
As a matter of fact, one of the most convincing evidence of saturation is the suppression of the away-side peak of di-hadron azimuthal correlations in p+A versus p+p collisions [8][9][10][11][12][13][14][15], due to the bigger target saturation momentum Q s in the p+A case. Recently, those calculations have been extended to the case of forward di-jets at the LHC, using an approximation of CGC expressions suitable for high-p T particles, dubbed improved transverse momentum dependent (TMD) factorization [16,17]. In this context, saturation effects are important when the total momentum of the di-jet system q T is of the order of Q s , and the first experimental measurements released recently [18] are compatible with that narrative.
In this kinematical regime however, the cross-section also receives important contributions from large logarithms ln 2 p 2 T q 2 T , that could reduce the impact of saturation effects. On the theoretical side, we know how to take into account the effects of the soft-gluon emissions responsible for those large logarithms. They affect the predictive power of the fixedorder perturbative expansion, and their resummation is therefore essential. The transverse momentum resummation formalism within the collinear factorization framework (C.F.) has been established for various processes, such as semiinclusive deep inelastic scatterings, Drell-Yan vector boson productions, and di-jet productions in hadronic collisions [19][20][21][22][23][24][25] in the past decades and has received great success in phenomenology. In particular, the transverse momentum resummation formalism gives good quantitative descriptions of transverse momentum distribution of heavy boson productions in the low-q T region at mid-rapidity [26][27][28][29][30][31][32][33][34][35][36][37]. Recently, it has been proven that such Sudakov resummation can be carried out together with the small-x resummation [38,39]. On the phenomenological side however, the interplay between the Sudakov and the small-x logarithms has only been studied in the context of low-p T di-hadrons [14], for which non-perturbative contributions dominate and prevent robust conclusions.
Before tackling the case of forward di-jets, we consider in this letter a simpler process, namely low-q T Z 0 -boson production, in order to study, for the first time, this interplay between Sudakov factors and non-linear small-x effects, and to assess to what extend saturation effects survive the soft-gluon emissions. We find that the suppression of the low-q T Z 0 -boson production in p+A versus p+p collisions is still present, but it gets significantly reduced. For completeness, we also compare our predictions with calculations obtained in the collinear factorization approach together with transverse momentum resummation, which can serve as a baseline without any small-x saturation effects. For that process, the largely-unknown collinear nuclear gluon PDFs at small-x play a key role, and as a result, the cross-section suffers from large uncertainties. Consequently, a measurement of the nuclear modification factor in agreement with our saturation prediction would also be compatible with collinear factorization with nuclear PDFs. Obviously, that conclusion could change drastically if precise small-x nuclear gluon distribution were to become available.
The plan of the letter is as follows. In section II we recall the calculation of the differential cross section for Z 0 -boson production in both the collinear factorization and CGC frameworks. In section III, we show our numerical results for the differential cross sections, compare with the available data, and give our predictions for the nuclear modification factor at forward rapidity. Section IV is devoted to conclusions.

II. THEORETICAL FRAMEWORKS
In our following calculations of Z 0 -boson production in proton-proton and proton-nucleus collisions, we adopt both the collinear and small-x framework together with the Sudakov type of transverse momentum resummation, and thus we can gain some quantitative insights into the onset of small-x effects through the comparison. These two frameworks are described in detail as follows.

A. Collinear factorization
In the transverse momentum resummation (also known as the Collins-Soper-Sterman resummation) formalism [21], the differential cross section of pp → Z 0 + X in the small transverse momentum region can be written as where the coupling is given by with θ W the Weinberg angle. The kinematic variables are defined as follows represents the collinear parton distribution function of parton a (b) inside the hadron A (B). For the purpose of numerical evaluations, we choose the parametrization of PDFs f a/A provided by CTEQ14 [40]. The Sudakov factor consists of two parts which read S sud (Q, b) = S p (Q, b) + S np (Q, b). At the next-to-leading logarithm (NLL) accuracy, the perturbative Sudakov factor S p (Q, b) is given by where the coefficients A 1 = C F αs(µ) 2π are derived from the one-loop correction, while A 2 = is computed from the two loop corrections of the double logarithmic term [41]. The QCD color factors are C F = 4/3, C A = 3 and T R = 1/2. Here we also follow the b * -prescription, This ensures that the integration over b ⊥ always stays in the perturbative region for the perturbative Sudakov factor. The non-perturbative Sudakov factor can be taken from BLNY parameterization [30,42,43] and it reads where the fitted parameters are found to be g 1 = 0.21 GeV 2 , g 2 = 0.68 GeV 2 , g 3 = −0.12 GeV 2 . In this transverse momentum resummation formalism, ) n , is the coefficient function to find parton i inside of parton a (i, j represent the flavor of quark. a, b could be either a quark or a gluon.), which can be calculated perturbatively [21,29,44].
i,g (z) = These coefficient functions take this simple form mainly because we have set the scale µ = µ b as in Ref. [29] and those terms proportional to ln µ/µ b vanish.
FIG. 1. Dominant diagrams contributing to the unintegrated quark/antiquark distribution function in the Drell-Yan process or Z 0 boson production in the CGC framework.
The dilute-dense framework uses the collinear PDF for the dilute projectile and the small-x gluon distributions for the dense target. The Sudakov resummation within the dilute-dense factorization framework at one-loop level has been demonstrated in [38,39]. It is the ideal framework to compute the Z 0 boson production in the forward rapidity at small-q T . In contrast with the di-hadron production in the forward rapidity [14], the non-perturbative effects are minimized, since fragmentation function is not involved in this process. The cross section in the dilute-dense factorization framework is given by Again the Sudakov factor S sud (Q, b) consists of two parts. The perturbative Sudakov factor is slightly different with that in the collinear factorization framework. The single logarithm term associate with the small-x target parton does not appear in the one-loop CGC calculation. Therefore only half of the B-term contributes to the perturbative Sudakov factor. Due to lack of the two-loop calculation, the corresponding A 2 -term in the dilute-dense factorization is not yet known. In this sense, our calculation in the dilute-dense factorization is only at the leading logarithm accuracy. Therefore, the perturbative Sudakov factor at one-loop level reads The non-perturbative Sudakov factor in the CGC framework, in principle, could be different with that in the collinear factorization framework. It can only be extracted from a global analysis to the experimental data. In this paper, we simply employ the same parameterization in the collinear factorization for simplicity.
, which can be computed from graphs as shown in Fig. 1, is the small-x unintegrated quark distribution in the coordinate space, which is similar with those derived in the DIS, semi-inclusive DIS and Drell-Yan processes [45][46][47][48][49][50][51][52][53][54] in the CGC framework. In coordinate space, this quark distribution reads is the modified Bessel function of the second kind. Phenomenologically, we use d 2 R ⊥ = 15 mb for the proton target and the rcBK solution [55] for the scattering amplitude N (x 2 , |b|). There are two α s 's in Eq. 8. The scale dependences can be different, although it is not quite obvious which scales should be used. In the phenomenological study, there are in general two prescriptions. (A) The α s in the numerator comes from the g → qq process and therefore its scale is chosen as µ b . However, the α s in the denominator comes from the definition of unintegrated gluon distribution function in the CGC framework [56], and its scale is related with the saturation scale. In our numerical calculation, we set the α s in the denominator to be 0.3. (B) Assume the scale dependences are the same, therefore they cancel each other. We find that prescription (B) fails to describe the experimental data. It is also worth noting that more extensive discussions on this issue can be found in Refs. [57,58]. We will return to this point in the next section.
In the experiments, the precision of the low-q T spectrum of the Z 0 -boson is limited by the energy resolution of the detectors [59][60][61]. The measurement with a finer bin-size is not achievable as well for the same reason. In order to make a better constraint on the theoretical inputs with the current experimental measurements, an optimized observable, φ * , which only involves the directions of momenta, is proposed in Refs. [59,60]. The resolution of the φ * measurement is supposed to be much better than the q T measurement, although the relevant physics is roughly the same. The rapidities of the leptons are denoted by η − and η + in the lab frame. After a Lorentz boost along the beam direction, we switch to a new frame where η − boost = −η + boost . As shown in Fig. 2, θ * is defined as the scattering angle between the lepton pair and the incoming protons in the new frame and ∆φ is the azimuthal angle of the di-lepton pair in the transverse plane. It is straightforward to get cos θ In the limit q T Q ≈ M Z , one finds ∆φ ∼ π and l 1 ≈ l 2 ≈ M Z 2 . Furthermore, we get sin θ * ≈ 2l T M and φ * = 2l T M Z tan π−∆φ 2 ≈ |q y T |/M Z , where q y T is the y component of q T as shown in Fig. 2. The φ * distribution is then given by, The case that φ * is small does not in general guarantee that q T is also small, since φ * is only determined by one component of q T . However, the cross section is dominated by small-q T region. Therefore the Sudakov resummation framework is still a good approximation. From Eqs. 1 and 6, one can also easily derive the corresponding expressions for φ * distribution under different frameworks by simply putting in a delta function according to the definition in Eq. 9.
At large φ * region, q T is forced to be large. First, we would need to implement the Y -term [21,29] to give the correct calculation for dσ/d 2 q T . Second, the power correction to φ * = M 2 ) becomes relevant and therefore Eq. 10 does not hold. In this paper, we only focus on the small-q T and small-φ * region.

III. NUMERICAL RESULTS
A. Z 0 boson production in the middle rapidity While the Z 0 boson is produced in the middle rapidity, the longitudinal momentum fractions of the incoming partons, x 1,2 ∼ M Z / √ S, are in general not small. This indicates that only the collinear factorization framework can be applied. The small-q T Z 0 boson cross section in the middle rapidity has been calculated numerically in several papers [26][27][28][29][30][31][32][33][34][35][36][37] under the collinear factorization and Sudakov resummation framework. In fact, the experimental data have also been included in the global analysis to extract the parameters in the non-perturbative Sudakov factor [30,42]. In this section, we present our numerical results in the middle rapidity and compare with the experimental data and previous studies in order to conduct a cross check.
Eq. 1 computes the Z 0 boson cross section, while, the experiments usually measure the di-lepton cross section at given rapidity window. Not all the Z 0 bosons produced in the middle rapidity will decay into di-lepton pairs that lie in the rapidity window covered by the detectors. In the Z 0 rest frame, the differential cross section of qq → Z 0 → l + l − is given by dσ/d cos θ ∝ A(1 + cos 2 θ) + B cos θ, with θ the angle between the outgoing lepton (l − ) and the incoming quark q. In general, both l − and l + are required to be at the same rapidity window in the experiments. The total di-lepton cross section is given by the integration of cos θ in the full phase space, which equals the Z 0 -boson cross section times decay fraction. The di-lepton cross section measured in the experiments is given by the cos θ-integral in the interval constrained by the rapidity window of the detectors. Therefore we get where x = cos[2 arctan e −η ] with η = min[y Z − y min , y max − y Z ]. y min and y max define the rapidity window of the di-lepton pair measured in experiments. 0.03366 is the fraction of the Z 0 boson decays to a certain flavor of di-lepton pair. When the rapidity window is wide, η is large and x ∼ 1. Most of the di-lepton pairs can be measured by the detectors. The cross section of di-lepton pair is simply that of the Z 0 -boson multiplies the decay fraction. However, in case that the rapidity window is very narrow, this factor could give a sizable suppression. qT distribution of di-lepton pair (through Z 0 boson) compared with data from CDF collaboration [62], D0 collaboration [63,64] and CMS collaboration [65].
We compute the differential cross section and compare with the experimental data measured by CDF [62], D0 [63,64] and CMS [65] collaborations in Fig. 3. It has already been shown in [30] the CDF data [62] is systematically 10% larger than the D0 data [63]. The parameterization for the non-perturbative Sudakov factor [30] used by us in the numerical evaluation is determined by fitting with the D0 data. Therefore, our results is about 10% smaller than the CDF data. Despite that, the theoretical results show good agreement with the experimental data. The shape of q T distribution is similar with that in Ref. [66], although different methods have been employed to deal with the non-perturbative effect. This confirms that the non-perturbative Sudakov factor has little impact on the shape of the q T distribution. But it may have some effects on the overall normalization. We also calculate the self-normalized φ * distribution in the collinear factorization framework using Eq. 10 and compare with the ATLAS data [61] in Fig. 4. As discussed in the previous section, the error band of the experimental data has been significantly reduced as compared to the q T measurement even with finer bin-sizes. The collinear framework calculation provides a good baseline description of the experimental data at small φ * .

B. Z 0 boson production in the forward rapidity
When the Z 0 boson is produced in the forward rapidity, the parton momentum fraction of the projectile side becomes large while that of the target side is small. At the LHC, √ S ≈ 5 TeV, x 2 is smaller than 10 −2 when y > 1. At Tevatron, √ S ≈ 2 TeV, x 2 < 10 −2 when y > 2. In these cases, the projectile proton is regarded as a dilute system while the target proton or nucleus can be viewed as a dense system. The so-called dilute-dense factorization becomes the ideal framework to study the Z 0 production process in the forward rapidity. On the other hand, the collinear factorization, which works perfectly in the middle rapidity, can still apply, since the momentum fraction of the target parton is not yet extremely small. In Fig. 5, we calculate the self-normalized Z 0 boson cross section in the forward pp collisions and compare with the D0 data [67]. In the CGC framework, there is no difference between proton and antiproton, since all the quarks and antiquarks come from the gluon splittings (i.e., they are sea partons.). The contribution from the valence anti-quark is missing. Therefore, the cross section calculated with the CGC framework is about 20% smaller than that given by the collinear factorization framework. However, this difference disappears in p+p collisions, as shown in the right two figures in Fig. 5.
As discussed in the previous section, the leftmost figure in Fig. 5 obviously shows that prescription (B) (the blue dotted curve) fails to yield the correct q T distribution, while prescription (A) along with the collinear factorization can. Therefore, it is essential to keep the different scale dependences in the couplings. In the paper, we will always adopt the prescription (A) in the numerical calculations.
The LHCb collaboration has measured the absolute differential cross section [68,69] in the forward p+p collisions as well. We compare our results with the experimental data in Fig. 6. Our theoretical calculation can describe the shapes of the distributions quite well, thanks to the success of the Sudakov resummation. However, we need to rescale our cross section by 1.15, which is larger than the one used in the middle rapidity, in order to match the absolute cross section. So far, we cannot understand well why a larger rescaling factor is required. However, according to Ref. [30], the overall normalization factor for different experiments in the BLNY parameterization could vary from 0.86 to 1.19. We assume this normalization factor remains the same in p+p and p+A collisions and therefore it is canceled in the nuclear modification factor. The modification of cross sections in forward p+p and p+A collisions is described by different languages in different frameworks. In the CGC framework, the gluon density is much larger in the nucleus than that in the proton and therefore the saturation effect is much stronger, thus R pA is smaller than one. In the collinear factorization approach, the cold nuclear effects are factorized out from the hard interaction and absorbed into the nuclear modified parton distribution functions (nPDFs). There are several parameterizations of nPDFs in the market [70][71][72][73][74][75]. In this paper, we employ the nPDFs provided by nCTEQ [73] and EPPS16 [75] groups to evaluate the nuclear modification factor, which is shown in Fig. 7, along with the prediction given by the CGC framework. In p+p collisions, the cross section  is self normalized in the range 0 < q T < 20 GeV, while in p+A collisions, it is normalized by 1 A 1 σpp dσpA dq T . By doing so, we can demonstrate not only the modification of the shapes, but also the suppression of the cross sections. The ratio of the "normalized" differential cross sections gives the nuclear modification factor, R pA = 1 A dσ pA /dq T dσpp/dq T . In collinear factorization, the error bands are given by uncertainties of nPDFs, while in the CGC framework, the error band is given by varying the saturation scale of the nucleus from 2Q 2 sp to 3Q 2 sp as the initial condition to solve the rcBK evolution equation. Fig. 7 shows that the CGC calculation agrees with the result given by the collinear factorization framework with EPPS16 nPDFs, albeit with a smaller error band. However, there is a clear difference between the CGC calculation and the nCTEQ15 nPDFs results. This is mainly because that the nPDFs at small-x are in general less constrained due to the lack of experimental data. The future measurements of the nuclear modification factor to the Z 0 boson production in the forward rapidity at the LHC should be included in the global analysis and would play an important role to reduce the uncertainty of the nPDFs at small-x. Furthermore, using Eq. 10, we can calculate the φ * distributions in forward p+p and p+A collisions which are shown in Fig. 8. In spite of a small difference between two frameworks in the calculation of dσ/dφ * , both CGC framework and collinear factorization framework with the EPPS16 nPDF predict the similar nuclear suppression factor, R pA . In addition, we also find that CTEQ15 nPDF predicts a stronger suppression at small-x.

IV. SUMMARY
To summarize, we calculated the differential Z 0 boson cross section in forward p+p and p+A collisions at small transverse momentum by applying Sudakov resummation in both collinear factorization framework and CGC framework. We find different nPDFs groups predict different nuclear modification effects and the uncertainties are very large at small-x. The CGC calculation coincides with the expectations given by EPPS16 nPDFs, with a relatively smaller error band. The upcoming LHC data can help to constrain the parameters in the nPDFs and improve our understanding about the dynamics of parton saturation, in particular the nonlinear evolution at small-x.