{'fraction_non_alphanumeric': 0.042780393977581745, 'fraction_numerical': 0.012513541597576883, 'mean_word_length': 4.035288878993654, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Wilson lines - color charge densities correlators and the production of eta' in the CGC for pp and pA collisions

We compute the inclusive differential cross section production of the pseudo-scalar meson eta' in high-energy proton-proton (pp) and proton-nucleus (pA) collisions. We use an effective coupling between gluons and eta' meson to derive a reduction formula that relates the eta' production to a field-strength tensor correlator. For pA collisions we take into account saturation effects on the nucleus side by using the Color Glass Condensate formalism to evaluate this correlator. We derive new results for Wilson line - color charges correlators in the McLerran-Venugopalan model needed in the computation of eta' production. The unintegrated parton distribution functions are used to characterize the gluon distribution inside protons. We show that in pp collisions, the cross section depends on the parametrization of unintegrated parton distribution functions and thus, it can be used to put constraints on these distributions. We also demonstrate that in pA collisions, the cross section is sensitive to saturation effects so it can be utilized to estimate the value of the saturation scale.

perturbative expansion. These have to be resummed and this leads to various factorization formalisms like the collinear factorization and the k ⊥ -factorization. In collinear factorization, meaningful physical quantities are obtained in terms of parton distribution functions (PDF). These distributions characterize the non-perturbative (large distance) physics and have to be determined experimentally from a fit to structure functions. This formalism can be applied on system where the typical exchanged momentum Q is hard, which means that it satisfies the inequality given by Λ 2 QCD ≪ Q 2 ∼ s where Λ QCD ≈ 0.2 GeV is the QCD scale and √ s is the center of mass energy. Note that Q 2 ∼ M 2 ⊥ ≡ M 2 + k 2 ⊥ , which means that unless the mass of the produced particle M is of the order of √ s, the collinear factorization is valid only for very large k ⊥ . This can be relaxed in the k ⊥ -factorization formalism which considers semihard collisions, meaning that the typical exchanged momentum obeys Λ 2 QCD ≪ Q 2 ≪ s. The resummation implemented in this approach takes care of large contributions that look like ln(Q 2 /Λ 2 QCD )α s n , ln(Q 2 /Λ 2 QCD ) ln(1/x)α s n and [ln(1/x)α s ] n [1,2,3,4]. This technique is used successfully to compute the production of many kinds of particles in high energy proton-proton collisions like heavy quarks [1,5,6,7,8] and in a number of other processes (see [9,10,11] for reviews of many applications). It is also used for predictions of Higgs boson production [12,13]. In this article, we use this formalism to compute the η ′ inclusive cross section at the RHIC energy.
When there is a nucleus involved in a collision at very high energy, there are new effects not included in the previous approaches due to the high density of gluons resulting from the emission enhancement at small-x (where x is the momentum fraction). These effects introduce a new scale Q s called the saturation scale at which the probability of having interactions between gluons of the same nucleus becomes important. At this transverse momentum scale, the gluons recombine and this slows down the growth of partons at smaller x. A naive estimation of Q s shows that it depends on x and the number of nucleons A like [14,15] so at small enough x or large enough A, the saturation scale is hard (Q 2 s ≫ Λ 2 QCD ). When the typical exchanged momentum is smaller than the saturation scale such as Λ 2 QCD ≪ Q 2 ≤ Q 2 s ≪ s, saturation effects have to be taken into account even if the system is still in the perturbative regime. This can be achieved in the Color Glass Condensate (CGC) formalism which is a semi-classical effective theory where the non- In this article, we are using the CGC to compute the inclusive differential cross section of η ′ meson in pA collisions at the RHIC energy ( √ s = 200 GeV). The main goal of this work is to look at the effect of saturation in η ′ production to validate the CGC approach and estimate the value of the saturation scale by comparing our predictions with experimental data. The η ′ is a pseudoscalar meson with a mass of M = 0.957 GeV, a decay width of 2 0.203 MeV and quantum numbers of I G (J P C ) = 0 + (0 −+ ) [16]. One of the most important features of η ′ is that it couples to the QCD anomaly [17,18]. One way to implement and model this physics is by introducing an effective interaction between gluons and η ′ mesons.
This was done in [19], where the authors are proposing a vertex that couples two gluons and a η ′ meson (ggη ′ ) to explain B mesons decay (B → η ′ + X). This vertex contains a form factor that depends generally on gluons and η ′ momenta and that can be related to the η ′ wavefunction. The structure of this vertex was investigated thoroughly using various techniques like the hard scattering and the running coupling approaches. [20,21,22,23,24,25]. We use these results on the gluons-η ′ coupling to study the η ′ production mechanism based on gluon fusion.
The η ′ production in pp collisions at high energy was studied in [26] where the exclusive cross section for the diffractive process p + p → p + p + η ′ is computed. In our study, we focus on the inclusive production mechanism p + p → η ′ + X which shares similar features with this previous analysis. The first attempt to compute η ′ production in high energy pA collisions was done by one of the present author in [27]. In this study, the collinear factorization is used to compute the cross section at RHIC by including intrinsic transverse momentum in the PDF with a Gaussian distribution. Based on physical arguments, the width of the Gaussian, which represents the typical transverse momentum of gluons inside the nucleus, is chosen to be Q s . The authors show that the η ′ production is sensitive to the saturation scale implemented in this way. However, they acknowledge that their calculation can be improved because they use the collinear formalism outside its range of validity. The goal of this article is to revisit the η ′ production with a more rigorous approach by doing a full CGC computation that includes recombination effects more realistically.
The computation of meson production in the CGC was undertaken in the past using mostly an "hybrid" approach where the proton and the nucleus are described by the collinear factorization and the CGC respectively. In this formalism based on pQCD-like techniques, the fragmentation function of collinear factorization is convoluted with the gluon or quark cross section computed in the CGC formalism. This is suitable for well-known mesons like pions for which a wealth of experimental data have been measured and for which fragmentation functions are well-known. Pion production for pA collisions is computed in [28,29,30] using this methodology. Contrary to pions, the data in high energy hadronic collisions for η ′ is scarce, so another approach is required to take care of hadronization effects and internal structure of the η ′ meson. In [31], an effective theory is used to estimate the tensor meson production in pp collisions. We use a similar approach in this article where the interaction between gluons and η ′ is described by an effective theory. As discussed previously, we include these effects in an effective vertex that includes a form factor. As will be shown in this article, this can be implemented easily in the CGC formalism.
We consider only the case of pp and pA collisions at RHIC. For nucleus-nucleus (AA) collisions, the total number of η ′ -mesons produced by semihard collisions in the first instants (for t < 1 fm/c) should be important and the saturation effects would also be present.
Experimentally however, they cannot be detected because most of them decay inside the medium created by the collision. This is because the η ′ mean lifetime, which is about t η ′ ≈ 4.93 fm/c, is smaller than the time where the medium exists, which is from 1 fm/c up to 10 fm/c. Moreover, by considering pp and pA, we avoid all the complications that would result from the creation of the medium which include the understanding and modelling of the quark-gluon plasma properties. Finally, there are analytical solutions for the gauge field in pp and pA collisions, while the analytical solution in AA is still elusive. For these reasons, our present analysis is only applied to pp and pA collisions.
This article is organized as follows. In section II we describe the effective vertex used throughout the rest of the article. In section III, we show how to compute η ′ production in the CGC for pA and pp collisions. We start by deriving a reduction formula that relates the cross section to a correlator of field-strength tensors. This correlator is then evaluated to leading order in pA collisions and the result can be interpreted in terms of physical processes.
We also show how the k ⊥ -factorized cross section for pp collisions can be recovered in the low density limit of pA cross section. In section IV, we compute the correlators appearing in the expression of the pA cross section using the McLerran-Venugopalan model. In section V, we evaluate numerically the cross section for pA and pp and discuss the range of validity of our computation. Sections III and IV contain a lot of technical details. The reader interested in results can jump directly to section V.
Throughout the article, we use both light-cone coordinates defined by and Minkowski coordinates. It should be clear by the context which one is used. We also use the metric convention g µν = (1, −1, −1, −1).

II. EFFECTIVE THEORY
The effective theory used in this article couples gluons and the η ′ meson. In momentum space, the g * g * η ′ effective vertex (where g * means off-shell gluon) is given by where ǫ µνρα is the Levi-Civitta antisymmetric tensor, M is the η ′ mass, a and b are color indices, p ρ and q α are gluon momenta and F (p 2 , q 2 , M 2 ) is the η ′ form factor. The explicit expression of the interaction vertex have been studied in a number of articles where different parametrizations of the form factor can be found [19,20,21,22,23,24,25]. To get a first approximation of η ′ production and because we are mostly interested in making a comparative study between pp and pA collisions, we use a simple expression given by [26] F (p 2 , q 2 , where H 0 = F (0, 0, M 2 ). To get a better approximations of η ′ production, other parametrizations should be used. In the limit of no gluon virtualities (p 2 , q 2 = 0), the form factor is a constant that can be fixed by looking at the decay of ψ → η ′ + γ. It is given by The main applications of this coupling are related mostly to B and Υ decay where processes such as g * → g + η ′ and g * + g * → η ′ are considered [19,32,33,34,35]. More recently, gluon fusion was used to compute η ′ production in high energy hadronic collisions [26,27,36] and from a thermalized medium [37].
It is convenient for our purpose to consider the interaction Lagrangian given by that reproduces the vertex Eq. (2) in the perturbative expansion. As seen in the next section, this can then be used to derive a reduction formula. Here, G µν a (x) is the usual field-strength tensor given by where A µ a is the gauge field of gluons andG µν a (x) = ǫ µνρσ G a,ρσ (x) is the dual field-strength tensor. The Lagrangian is non-local because the vertex includes a form factor. It can be easily seen that L int (x) contains three types of vertices, namely ggη ′ , gggη ′ and ggggη ′ . At leading order however, only the first one is necessary and considered in this article.

III. PRODUCTION OF η ′ FROM THE CGC
In collisions at very high energy, the wave function of nuclei is dominated by gluons that have small longitudinal momenta (soft gluons) because of the emission enhancement at small-x. The CGC is a semi-classical formalism that describes the dynamics of these degrees of freedom. In this approach, the hard partons, which carry most of the longitudinal momentum, and soft gluons which have small longitudinal components, are treated differently. Because the occupation number of the soft gluons is large, classical field equations can be employed to understand their dynamics. The hard partons act as sources for these 5 classical field and are no longer interacting with the rest of the system (for reviews of CGC, see [14,15,38]).
In this formalism, computing a physical quantity involves two main steps. The first one is to solve the Yang-Mills equation of motion where the current J ν a (x) = δ ν+ δ(x − )ρ p,a (x ⊥ ) + δ ν− δ(x + )ρ A,a (x ⊥ ) represents random static sources localized on the light-cone [14,15] and D µ = ∂ µ − igA µ is the covariant derivative.
The functions ρ p,A (x ⊥ ) are color charge densities in the transverse plane of the proton and nucleus respectively. The next step is to take the average over the distribution of color charge densities in the nuclei with weight functionals W p,A [ρ p,A ]. For any operator that can be related to color charge densities, this can be written as Computing the weight functional is a highly non-perturbative procedure so it usually involves approximations based on physical modelling. In the limit of a large nuclei at not too small x, it can be approximated by the McLerran-Venugopalan (MV) model, which assumes that the partons are independent sources of color charge [39,40]. Within this assumption, the weight is a x A independent Gaussian distribution and the two point correlator is simply [14,15,39,40] where µ 2 A = A/2πR 2 is the average color charge density and R is the radius of the nucleus. It is assumed here that the nucleus has an infinite transverse extent with a constant charge distribution. Edge effects can be included by changing µ 2 → µ 2 (x ⊥ ) and by choosing a suitable transverse profile. Throughout this article, we only consider the constant distribution case.
Within the MV model, the weight functional does not depend on longitudinal coordinates and therefore, the model is boost invariant. This however can be relaxed by considering the quantum version of the CGC. In that theory, quantum radiative corrections become important below a certain scale x 0 ≈ 0.01. These corrections can be resummed by using a renormalization group technique which leads to the JIMWLK evolution equation [41,42,43,44,45]. In the quantum CGC, the weight functionals W 1,2 [ρ 1,2 ] obey this non-linear evolution equation in x. Because the MV model is valid in the range x ≈ 0.01 − 0.1, it can be used as an initial condition for the evolution at smaller x. In this article however, we consider only the regime where the MV model is valid and do not consider the small-x evolution although it could be done in principle. 6 On the proton side, the average computed with W p [x p , ρ p ] can be related to the unintegrated parton distribution function (uPDF) φ 1 like where N c is the number of color. By construction, the uPDF obeys and is normalized such that where xG(x, µ 2 ) is the collinear parton distribution function and µ 2 is the factorization scale.
The uPDF can be obtained from a fit to structure function and evolved to the desired value of x p , Q 2 and p 2 ⊥ using evolution equations such as the BFKL or the CCFM equations. One important ingredient is missing for the computation of η ′ meson production cross section. We need a relation between the cross section and a correlator that can be evaluated using the CGC formalism. This is done in the next section using a reduction formula and the effective theory.

A. Reduction Formula and the Cross Section
The computation of η ′ mesons from the CGC can be calculated from a reduction formula.
The starting point is the expression of the average number of η ′ produced per collisions given byn = ∞ n=1 nP n where P n is the probability to produce n particles. This can be converted to an equation in terms of creation/annihilation operators that can be evaluated in quantum field theory. This is given by [46] where |0 in is the in vacuum. Then, the standard LSZ procedure can be used to write this as [47] (2π where Z is the wave function normalization and where we assumed the asymptotic conditions lim t→±∞η (x) = √ Zη out,in (x) for the η ′ field operator in Heisenberg representation. It is possible to use the equation of motion ofη(x) given simply by 7 where we definedT (x) ≡ 1 Ĝ µν a (y)ˆ G a,µν (z) to rewrite the reduction formula in a convenient way. We get finally that whereT (k) is the Fourier transform ofT (x) evaluated at a η ′ meson on-shell momentum and where we set Z = 1. The angular brackets Ô here indicates expectation value ofÔ in the initial state.
The only assumptions used in deriving Eq. (15) are: • There are no η ′ mesons in the initial state.
• The η ′ is produced on-shell.
The first assumption is justified by the fact that in high-energy collisions, the number of η ′ in a hadron before the collision (in the initial state) is negligible. This allows us to use the in vacuum and the fact that a in (k)|0 in = 0 to simplify the reduction formula.
Using the second assumption, we can treat the η ′ meson as a stable particle which can be produced on-shell and which is well described by the free spectral density that looks . Therefore, by making this assumption, it is possible to use the asymptotic conditions described earlier. However, η ′ -mesons are resonances, so the spectral density should look rather as a Breit-Wigner function ρ( where Γ is the decay width. These effects however are taken into account by the form factor F (q, p, M).
Eq. (15) is the main result of this section. It relates the average number of η ′ mesons produced to a correlator of field strength tensors. This correlator can then be evaluated using any analytical or numerical methods. The average ... depends on the system studied.
Looking at a plasma at equilibrium, it could be computed using finite temperature field theory or the AdS/CFT correspondence. These two formalisms are relevant to nucleusnucleus collisions where a medium at equilibrium is created. We are interested here in pA and pp collisions where no such medium is formed so these techniques are not pursued in this study. Rather, we use the CGC which describes initial state and saturation effects in high-energy hadronic collisions.
Having expressed the average number of η ′ produced in terms of a correlator of field strength tensors, it is possible to compute the inclusive cross section in the CGC formalism 8 which is given by [46,48] In this expression, b ⊥ is the impact parameter. The fields T are functionals of the source once the Yang-Mills equation of motion of the gauge field is solved (see Eq. (5) for the expression of the field strength tensor as a function of the gauge field).

B. Cross section in pA Collisions
In pA collisions, there are two saturation scales (one for the proton (Q p ) and one for the nucleus (Q A )) that satisfy Q p < Q A . When the transverse momentum of the η ′ is small enough, the nucleus is in a saturation state while the proton is not because we have In that case, the system is semi-dilute, meaning that one of the source is strong (or equivalently, the typical transverse momentum is small) and obeys ρ A,a (k ⊥ )/k 2 ⊥ ∼ 1 while the other source is still weak ρ p,a (k ⊥ )/k 2 ⊥ ≪ 1 [38,49]. The weak source can be used as a small parameter to solve the Yang-Mills equation perturbatively. The solution of the gauge field can be computed analytically to all orders in ρ A,a (k ⊥ )/k 2 ⊥ and to first order in ρ p,a (k ⊥ )/k 2 ⊥ in different gauges [50,51,52,53]. We use here the solution in the light-cone gauge of the proton [52] but in Appendix C, we perform the same calculation in covariant gauge to show that our result is gauge invariant.

Gauge Field and Power Counting
In the light-cone gauge of the proton (A + = 0) with a nucleus in covariant gauge moving in the negative z direction, the solution of the gauge field in pA collisions can be separated [52]. Note that the field The explicit solution is given by [52,54] A i p,a (k) = 2πgδ(k − ) where U ab (k ⊥ ) is a Wilson line in adjoint representation defined in Eq. (35), g is the usual QCD coupling constant, δ ab is the Kronecker delta in color space and f abc is the structure is non-zero and is related to A i pA,a (k) but it does not appear in the final expression of the cross section so it is not needed in the computation of η ′ production. All the other components are zero.
The production cross section of η ′ mesons is related to a field strength tensor correlator given by where p 2 , q 2 = k − p, q. In principle, a correlator like this contains contributions from all orders in both sources. Because the proton source is weak, it is possible to simplify this considerably using a power counting argument to isolate the leading order contribution.
For the sake of this power counting argument, we use a n (...) which denotes terms having n gauge fields A µ (...) and where a A ∼ O(ρ A ) and a p , a pA ∼ O(ρ p ). At first, let us consider only the contributions from the abelian part of the field-strength tensor. The terms in these contributions have four powers of gauge field such as B abelian ∼ a 4 tot ∼ (a A + a p + a pA ) 4 . Naively, one would expect the dominant contribution to come from terms that have many powers of the nucleus gauge field like a 4 However, these terms vanish because of the Lorentz structure. For example, a typical term in the abelian contribution would look like T abelian ∼ ǫ µνρσ p µ q ν A ρ A σ . When we sum on indices, this kind of term will contain at most one strong gauge field A − ∼ a A . Thus, the dominant . Using a similar argument, it is possible to show that the non-abelian part have no leading order contribution in the sense that it is for example also vanish because of the Lorentz structure of the correlator. This is because the typical non-abelian contributions look like Once the Lorentz indices are summed, the second typical term T ′′ = 0 because in this gauge, A + = 0. For T ′ , it contains only one strong field A − ∼ a pA but it contains two powers of weak field like A i ∼ a p , a pA . Thus, when it is squared, it gives at most a contribution of B non−abelian ∼ O(ρ 2 A ρ 3 p ).

C. Evaluation of the Correlator
Using the explicit expression of the field strength tensor in terms of gauge field and keeping only the dominant and non-zero contributions, the correlator can be written as where ǫ ij is the Levi-Civitta antisymmetric tensor with i, j = 1, 2 and where we defined To obtain the preceding expression we make the assumption that the virtualities in the form factors are due solely to transverse momentum such as . This approximation is necessary to recover k ⊥ -factorization in the dilute limit as seen in section III E.
It is convenient to separate B(k) in four different terms such as where z, z ′ = {p, pA} and where we performed the integration on p + and q + using the delta functions in Eq. (18). These four terms can be evaluated explicitly by substituting the solution of gauge fields Eqs. (17), (18) and (19).
For the first term, it is a straightforward calculation to show that The calculation of the second term is similar but requires some more work. By direct substitution of the expression of the gauge fields, we have The integration on the longitudinal momentum p − can be done by looking at the analytical structure of the equation. First, we write a part of the integrand as In the complex plane of p − , the first term in the RHS of Eq. (25) has two poles on the same side of the real axis and goes like 1 (p − ) 2 when p − → ∞. Thus, closing the integration contour in the upper-half plane and using the residue theorem, the integration on p − of this term leads to a zero contribution because the contour at infinity has a zero contribution and because the contour does not enclose any singularities. For the second term of Eq. (25), we use the principal part (PP) identity 1 x±iǫ = PP 1 x ∓ iπδ(x). The integration on the principal part is zero because the integrand does not depend on p − and lim R→∞ R −R dpPP 1 p−a = 0 while the delta function integration is trivial. We finally get The calculations of B pA,p (k) and B pA,pA (k) are similar. Going through the same steps as for and The final result for the correlator B(k) = B p,p (k) + B p,pA (k) + B pA,p (k) + B pA,pA (k) can be written compactly as where we definedŨ From this equation for B(k), we can now evaluate the differential cross section. It is given Thus, we can relate the η ′ production cross section to a correlator of sources and Wilson lines. This will be evaluated in the next section in the MV model. Clearly, this expression does not have the k ⊥ -factorization structure but it can be recovered in the dilute limit. First, we can look at the physical interpretation of the different terms in the cross section.

D. Diagrammatic Content and Physical Interpretation
To facilitate the physical interpretation, it is convenient to interpret the gauge field expressed in Eqs. (17), (18) and (19) in terms of Feynman diagrams as shown in Fig. 1.
The strong gauge field A µ A corresponds to a resummation of tree diagrams with any number of strong source insertions [55]. The weak gauge fields A µ p and A µ pA also resum an infinite number of tree diagrams, the difference being that they contain one weak source insertion [50]. Then, the η ′ production and the correlator given in Eq. (21) can be represented diagrammatically in Fig. 2 and 3. These figures show all the diagrams included in the calculation.
The first term in the figure corresponds to the part where the η ′ is produced from gluon fusion. The second term contains multiscattering effects and eventually, saturation effects.
Overall, this leads to the following physical picture. A gluon inside the proton interact with the classical background field of the nucleus and gets multiscattered. Once it has gone through the nucleus, it combines with a gluon and produce the η ′ .

E. Recovering k ⊥ -factorization in the cross section
It is possible to recover k ⊥ -factorization from Eq. (31) by looking at the dilute limit of the nucleus characterized by a weak source such as ρ A,a ≪ 1. In that case, we are allowed to keep only the first term of the Wilson line expansionŨ ae (k ⊥ ) = 2(2π) 2 δ 2 (k ⊥ )δ ae + O(ρ A ).
In this low-density limit, the cross section becomes In terms of Feynman diagrams, this expression corresponds to neglecting all the multiscattering diagrams shown in Fig. 3. The neglected diagrams are the ones that break k ⊥ -factorization as can be seen from the following argument.
In this dilute limit, the correlator of the nucleus can also be related to the uPDF φ 2 like in Eq. (9). Using this results, we find that the cross section is given by This is the k ⊥ -factorized expression of the cross section and is totally equivalent to Eq.
(A7) obtained directly from the k ⊥ -factorization formalism (see Appendix A). Thus, in the low-density limit of pA collisions, we recover a formalism that describes pp collisions in the semihard regime. This is very similar to quark and gluon production in pA collisions [48,50,56,57]. Note that to obtain this result, it is necessary to assume from the beginning that the form factors depend only on transverse momenta.

IV. COMPUTATION OF CORRELATION FUNCTIONS IN THE MV MODEL
In this section, we compute the relevant correlators appearing in our expression of the cross section using the McLerran-Venugopalan (MV) model. Throughout this calculation, we use the notation of [58,59]. We are interested in correlators containing both Wilson lines and color charge densities such as the ones included in Eq. (31). In App. D, we discuss the general case and give more details on the calculation. Note here that to make sense of the ordered path, we start with color charge densities that depend on the longitudinal coordinate x + . At the end of the calculation, we take ρ a (x + , x ⊥ ) = δ(x + )ρ a (x ⊥ ) since we use the MV model which assumes that the nucleus is moving at the speed of light. In this more general case, the 2-point function is simply where µ 2 (x + ) is the average color charge density at point x + . It is related to the average color charge density by µ 2 = dx + µ 2 (x + ) = A/2πR 2 where R is the radius of the nucleus.
In this model, W [ρ] is still Gaussian, so all even-point functions can be written in terms of the 2-point function using Wick theorem and all odd-point functions are zero.
The Wilson line is defined as where t c are the SU(N c ) generators in adjoint representation, P + is the path ordering in the light-cone coordinate z + and G 0 is a Green function solution of With these definitions, it is possible to compute the needed correlators. The general strategy is to express the correlators in terms of the following known results for Wilson lines in adjoint representation [59]: and where we defined the quantityμ 2 (b + , a + ) ≡ b + a + dz + µ 2 (z + ). In the next subsections, we express the correlators appearing in the η ′ cross section in terms ofŪ andV which are defined in Eqs. (37) and (38).

A. 1 Wilson line -1 color charge density correlator: the main building block
The first correlator does not appear explicitly in the η ′ production cross section but it is useful to understand the other results. We want to evaluate where we assume that b + > y + 1 > a + . For clarity, we define the sources included in Wilson lines as internal sources as opposed to external sources which appear explicitly in the correlator (like ρ c 1 (y + 1 , y 1⊥ ) in Eq. (40)). The first type is a contraction like ρ a j−1 (z + j−1 , z j−1⊥ )ρ a j+1 (z + j+1 , z j+1⊥ ) ρ c (y + , y ⊥ )ρ a j (z + j , z j⊥ ) where ρ c is an external source. The other ones have only internal sources. The second . Only the last one has a support and thus, a non-zero contribution.
The first step is to expand the Wilson line. The expression can then be written as where we used t a bc = −if abc and where f abc is the antisymmetric SU(N c ) structure constant. Using Wick theorem, the n-point correlation function can be expressed in terms of 2-point functions as the sum of all possible contractions. At first, we look only at the contractions of ρ c 1 (y + 1 , y 1⊥ ) with ρ a 1 (z + 1 , z 1⊥ )...ρ an (z + n , z n⊥ ). Using the MV expression for the two point function Eq. (34), we get As argued in [58,59], only adjacent sources can be contracted due to the ordering in z + . All the other nested and overlapping contractions have no support. It can also be shown that the contraction ρ a j−1 (z + j−1 , z j−1⊥ )ρ a j+1 (z + j+1 , z j+1⊥ ) (where ρ a j (z + j , z j⊥ ) is contracted with the external source) does not have support either. These properties are shown diagrammatically in Fig. 4. They can be used to split the correlator in two parts like By using these properties and by reorganizing the series, we have This complicated expression is just a combination of Wilson lines that is given more succinctly as This can be simplified further by using Eq. (37) and the fact thatŪ is an exponential. We can get easily that Specializing to the case of a charge distribution moving at the speed of light we have that ρ a (x + , x ⊥ ) = δ(x + )ρ a (x ⊥ ). Integrating both sides by y + , we get

B. 1 Wilson line -2 color charges densities correlator
In this subsection, we compute the correlator It is possible to devise diagrammatic rules shown in Fig. 5 that can be used to write F 2,1 (b + , a + ) in terms of known quantities. The rules are discussed in more details and  [58,59]. Finally, in this notation, the light-cone coordinates are ordered following the arrow, from left (the smallest) to right (the biggest) so that b + > a + .  generalized to all cases in App. D. Using these results, F 2,1 (b + , a + ) can be represented diagrammatically as in Fig. 6. This includes all possible contractions and topologies that need to be resummed. According to the rules of App. D, it can be written as This can be simplified further by using the explicit expressions given in Eqs. (34), (37) and (46). We get that Considering that the nuclei is moving at the speed of light and integrating on both sides by y + 1 and y + 2 , we get where we defined

C. 2 Wilson lines -2 color charges densities correlator
In this subsection, we compute the correlator This calculation is similar to the one of the second correlator F 2,1 using the diagrammatic rules. Three typical diagrams of F 2,2 (b + , a + ) out of seven are shown in Fig. 7. The first one corresponds to the contraction of the two external sources and the other ones corresponds the connected part which consists in all contractions between internal and external sources.
The first diagram is straightforward to compute. It is given by Using Eqs. (34) and (38), we get that The second diagram shown in Fig. 7 is given by Using Eqs. (37) and (38), we get which is zero because of the color structure. All the other diagrams included in F 2,2 (b + , a + ) can be computed in a similar way. There are five other different ways of inserting the source and they all vanish because of the color structure. The only non-zero term is the first one and so we have F 2,2 (b + , a + ) = F 2,2 1 (b + , a + ). Finally, with a nucleus moving at the speed of light, we obtain This concludes the computation of correlators. We are now in a position to evaluate the cross section within the MV model. To make a comparison with pA collisions for a wide range of transverse momentum where the MV model is valid (see Fig. (9)), we chose the rapidity y = 1. The results of the numerical calculation are shown in Fig. (8) where we also present the result for the pA case at small saturation scale (Q 2 s = 1 GeV) and with the proton described by the parametrization CCFM J2003 set 3.
Then, we can write To obtain this expression, we used Eq. (9) to convert the average on proton sources to an unintegrated distribution function and we Fourier transformed the correlator in the second term. We also neglected all terms with only one Wilson line because they are numerically very small compared to the other terms. This can be seen as follows. First, using translation invariance of the correlator, we see that averages with one Wilson line in the cross section are proportional toŪ(0) (see Eq. (62)). This quantity is small because there is an infrared singularity appearing in the argument of the exponential in the following way [59]: This is not exactly zero because the infrared singularity is regulated by non-perturbative effects (confinement) and these induce a cutoff of order Λ QCD . The saturation scale in a large nuclei at small-x satisfies Q 2 Thus, the correlators with one Wilson line can be neglected.
It is then a straightforward calculation, involving some change of variables, translation invariance and Eqs. (34) and (61), to obtain where S ⊥ = (2π) 2 δ 2 (0) is interpreted as the transverse area of the nuclei. We also introduced here an infrared regulator given by Λ in one of the denominator 1 . This is because gluons separated by a distance greater than the size of a nucleon (∼ 1 fm) are not correlated because of confinement so that Λ = Λ QCD . This is not taken into account explicitly in the MV model, leading to infrared divergent quantities [67] which are regulated by adding the infrared regulator [56]. This is equivalent to define a gluon distribution function for the MV model as We do not introduce a regulator in the other denominator 1 p 2 ⊥ because as seen in Eq. (A7), it is part of the matrix element |M| 2 ∼ ǫ ij ǫ kl which is well behaved in the infrared.
In the first term of the cross section, one integration can be done easily with the delta function. The remaining integrals will be done numerically. The second term can be simplified further by doing some integrals analytically as shown in Appendix B. The first term is very similar to the pp cross section and does not involve any saturation scale dependence because it does not include any multiscatterings. The Q s dependence occurs in the second term through the expression ofV where we define it as . (67)

Kinematical range
The cross section derived in the previous section is restricted to a certain range of validity because we are using the MV model. This model can be used when the radiative corrections which goes like α s ln(1/x) are not too large, which is when 0.01 production, the momentum fraction of gluons in the proton and the nucleus are given by The kinematic range in terms of the η ′ transverse momentum and rapidity where 0.01 x A 0.1 is depicted in Fig. 9 for RHIC energies.
It is possible to extend the range of validity to smaller values of x A by using the JIMWLK equation [41,42,43,44,45]. This renormalization group equation resums the large radiative corrections due to very small x modes. The effect is to change the correlation between sources and in general, one looses the Gaussian structure of the weight functional W [ρ].
This complicates the computation of Wilson line correlators and is outside the scope of this article.

Numerical Results
We present in this section the numerical results for η ′ production. We integrate numeri- In Fig. (12), we show the result for the inverse nuclear modification factor R pA for y = 1 defined as We use this because the pp cross section is zero at |k ⊥ | = 0. for RHIC energies in terms of the η ′ transverse momentum k = |k ⊥ | and rapidity y.

C. Analysis
For the pp cross section, there is a wide range of variability for the cross section depending on the parametrization of the uPDF used. This feature can be used to discriminate between the different parametrizations to determine the most accurate one. Thus, by combining this analysis with experimental data, the η ′ production becomes another observable that can be utilized to constrain models of uPDF. A similar conclusion was reached in [31] where the f 2 -meson production was studied. Note however that the results of the cross section at small transverse momentum should be treated carefully. As stated in the appendix A, the validity of the k ⊥ -factorization approach depends on the presence of a large scale compared to the QCD scale given in our case by Q 2 ∼ M 2 ⊥ . At very small momentum, we have that the following way. The saturation scale is defined as the momentum at which the probability of interaction between different parton cascades is of order one [14,15]. The partons having a transverse size satisfying δx ⊥ ∼ 1 Q ⊥ ≥ 1 Qs will have a very high probability to recombine. As the saturation scale is increased, there will be more gluons that will be sensitive to these nonlinear effects. The final result is that as more and more gluons have a high probability to recombine, the growth of their population will decrease and so is the cross section for a given set of parameters. This effect is clearly seen in the cross section computed in our model.
The fact that the η ′ production cross section is sensitive to the saturation scale can be used to estimate the numerical value of Q s . The current estimate based on HERA data for deep inelastic scattering shows that for RHIC energy, it is given by Q s ∼ 1 − 2 GeV [14,15]. It would be interesting to compare our calculation to experimental data to make an estimation of Q s based on η ′ production and see if it is consistent with the previous estimation.
The nuclear modification factor shows clearly the saturation effects for 0.4 |k ⊥ | 3.4 because in that range, R pA > 1. At large transverse momentum, it approaches the value R pA ≈ 0.45. One would expect R pA to be one in that range of momentum because it corresponds to the regime where there are no saturation effects and where the k ⊥ -factorized cross section can be used. This discrepancy can be explained by looking at the approximations made during the calculation. We neglect two terms in the cross section that are proportional to one Wilson line (see Eq. (64)). In the dilute limit (when |k ⊥ | > Q s ), these terms cannot be neglected and have a non zero contribution that would make R pA ≈ 1. In that sense, our approximation fails at large transverse momentum and our result should be treated carefully in that regime.

VI. CONCLUSION
In this article, the inclusive cross section for η ′ production in pp and pA collisions is computed. The pp case is analyzed for different reasons. First, it serves as a basis to validate our model for η ′ production against experimental data. This model includes two main ingredients. The first one is the effective theory which is at the foundation of our analysis. We use in our study a very simplified form of this effective theory where the vertex is given by Eq. (3). As shown in [19,20,21,22,23,24,25], this is a very crude approximation of the real vertex, so there are some improvements that can be done in this direction in the future. The second ingredient are the uPDF parametrizations. As seen in Fig. (8), there is still a large variability in the predictions made by different uPDF.
Therefore, η ′ production could be used to constrain the models of uPDF once it is compared to experimental data. The same conclusion was obtained in [26] for the exclusive process p+p → p+p+η ′ . Finally, from the theoretical point of view, we studied pp collisions to see if the cross section can be obtained as the low density limit of the pA cross section. There are now many known examples where this can be seen like gluon production [50,56,68], quark production [48,57] and tensor meson production [31]. We have shown that η ′ production also obeys this property and in that sense, it is a consistency check for the approach used for the pA case.
The pA results for the η ′ inclusive cross section show some very interesting features. First, we show that they are sensitive to the value of the saturation scale. This can also be seen in the plot of the nuclear modification factor. This property could be used to make an estimate of Q s by comparing with experimental data, which is one of the main goal of our analysis.
This information is very important for the study of other particle production such as quark and gluon production in both pA and AA where saturation effects play an important role.
Thus, a measurement of η ′ at RHIC would improve our knowledge of gluon distribution in a nucleus.
Throughout the article, we assumed that the gluon fusion process was the dominant mechanism in η ′ production. There is one other production process that could also be important. It is the photon fusion where two off-shell photons emitted by the protons or the 30 nucleus merge to give a η ′ such as γ * + γ * → η ′ + X. This is estimated in [26] and according to this analysis, it should be subdominant in the exclusive production. We assume that this holds also in our study, although a careful analysis of this process should be performed. Of course, our calculation could be made more accurate by investigating this last issue. where we assume for simplicity that b + > y + 1 > y + 2 > ... > y + m−1 > y + m > a + . To get the general result where all the y + j can have any value between b + and a + , one has to sum over all orderings. All other cases where a source is missing such as in the definition Eq. (D4) can be treated in a similar way. It is now convenient to define a new set of variables c + j which have the following property y + j > c + j > y + j+1 . Then, we use these new variables to split the + coordinate in slices, which allow us to write the Wilson lines (using their properties) as U a j b j (b + , a + |x j⊥ ) = U a j α (j) (c + p , c + p+1 |x j⊥ ) U α (j) m+1 b j (c + m−1 , a + |x j⊥ ). The only unknown quantity left is H 1,n (b + , a + |{a}, {b}), to which we now turn. Remember that in this quantity, there is only one external source and this source, using Wick theorem, can be contracted with an internal source in any of the n Wilson line included in H 1,n (b + , a + |{a}, {b}). In Sec. IV A, we show that when a source is contracted with a Wilson line, it breaks the Wilson line in two, change the color structure and multiply the overall expression by µ 2 (y + )G 0 (x ⊥ ). For the general case, it is shown in Fig. (14). Using this and the definition of H, we get H 1,n (b + , c + 1 |{a}, {b}) = n j=1 µ 2 (y + )G 0 (y ⊥ − x j⊥ )f cdd ′ × F n (b + , y + |{a}, {β}) β j =d F n (y + , a + |{β}, {b}) β j =d ′ . (D10) which completes our calculation since the F n 's are assumed to be known from the beginning.
This general result can be summarized into diagrammatic rules to compute any F m,n .
• Determine a +-ordering for the external sources.
• Draw all possible ways the external sources can be contracted with the Wilson lines and with each other by respecting the +-ordering.
• A blob with n external line is given by F j (b + , a + |{a}, {b}).
• A source insertion at y + with a color index c on line j multiplies the expression by µ 2 (y + )G 0 (y ⊥ − x j⊥ )f cdd ′ and changes the indices on the left blobs to β (j) = d and on the right blob to β (j) = d ′ .