Linearly-polarized small-x gluons in forward heavy-quark pair production

We use the Color Glass Condensate (CGC) framework to study the production of forward heavy quark-antiquark pairs in unpolarized proton-nucleus or proton-proton collisions in the small-x regime. In the limit of nearly back-to-back jets, the CGC result simplifies into the transverse-momentum dependent (TMD) factorization approach. For massless quarks, the TMD factorization formula involves three unpolarized gluon TMDs: the Weizs\"{a}cker-Williams gluon distribution, the adjoint-dipole gluon distribution, and an additional one. When quark masses are kept non-zero, three new gluon TMDs appear -- each partnered to one of the aforementioned distributions -- which describe the distribution of linearly-polarized gluons in the unpolarized small-x target. We show how these six gluon TMDs emerge from the CGC formulation and we determine their expressions in terms of Wilson line correlators. We calculate them analytically in the McLerran-Venugopalan model, and further evolve them towards smaller values of x using a numerical implementation of JIMWLK evolution.


I. INTRODUCTION
In hadronic reactions that are governed by more than one hard momentum scale, the standard QCD framework of collinear factorization at leading twist becomes insufficient, and one needs to resort to more sophisticated factorization schemes. One such scheme is TMD factorization [1][2][3][4][5][6][7][8], which makes use of transverse-momentum-dependent parton distributions, or TMDs for short. One of the many intricacies of TMDs is the fact that, in contrast to the usual collinear PDFs, their operator definition depends on the hard process under consideration, hence at first glance, universality is broken.
More specifically, forward quark-antiquark pair production in dilute-dense collisions is characterized by three momentum scales: P t , the typical transverse momentum of a single quark, and always one of the largest scales; k t , the total transverse momentum of the pair, which is a measure of the transverse momentum of the small-x gluons coming from the target; and Q s , the saturation scale of the nucleus, which is always one of the softest scales. The value of k t with respect to Q s and P t governs which factorization scheme is relevant. Indeed, when k t ∼ Q s P t (the quark and the antiquark are almost back-to-back), there are effectively two strongly ordered scales k t and P t in the problem and TMD factorization applies [12], implying the involvement of several gluon TMDs that differ significantly from each other, especially in the saturation regime, when k t ≤ Q s [22]. In the other regime: Q s k t ∼ P t , k t and P t are of the same order and far above the saturation scale, hence high-energy factorization [34,35] is applicable. In this case, only the linear small-x dynamics governed by the Balitsky-Fadin-Kuraev-Lipatov (BFKL) equation [36][37][38] is important, and the TMDs differ no more, implying that only one such distribution plays a role. Interestingly, both regimes, i.e. the TMD regime and the high-energy factorization regime, are encompassed within the CGC approach [16].
Indeed, in [12,39], the cross section for forward di-jet production in proton-nucleus collisions was calculated within the CGC. It was then shown that, in the back-to-back limit k t ∼ Q s P t , a TMD factorization formula could be extracted, the result being the same as in a direct TMD approach (i.e., without resorting to the CGC). However, in contrast to the direct TMD approach, the calculation in the CGC yields explicit expressions for the TMDs in terms of Wilson lines, which can be evolved in rapidity through the nonlinear Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) equation, as was demonstrated in [22].
In this paper, we build further on that work, by studying the forward production of a heavy quark-antiquark pair. As already observed earlier (see for instance [2,11,14,15,40,41]), by keeping a non-zero quark mass, the cross section becomes sensitive to additional TMDs, which describe the linearly-polarized gluon content of the unpolarized target, or in our case, nucleus. The three unpolarized gluon TMDs that describe the gluon channel gA → qq will be accompanied by three 'polarized' partners, which couple through the quark mass and via a cos(2φ) modulation, where φ relates to the quark-antiquark pair and is defined below. This is analogous to what happens in the γ * A → qq process (in that case not only a non-zero quark mass but also a non-zero photon virtuality brings sensitivity to linearly-polarized gluons), although there only one unpolarized gluon TMDs is involved (the Weizsäcker-Williams distribution), along with its polarized partner [11,14].
The paper is organized as follows. In section II, we give the result of the CGC calculation of the forward heavy quark-antiquark pair production cross section, and demonstrate how the six gluon TMDs (three unpolarized and three linearly-polarized) emerge in the appropriate limit. In section III, we compute these TMDs analytically in the McLerran-Venugopalan (MV) model and compare our results with the existing literature, after which in section IV they are numerically evaluated and evolved in rapidity with the help of a lattice implementation of the JIMWLK equation. Finally, we conclude and give an outlook for further work.

II. EXTRACTING A TMD FACTORIZATION FORMULA FROM THE CGC FRAMEWORK
We consider inclusive quark-antiquark pair production in the forward region, in collisions of dilute and dense systems The four-momenta of the projectile and the target are massless and purely longitudinal. In terms of the light-cone variables, x ± = (x 0 ± x 3 )/ √ 2, they take the simple form p p = s/2 (1, 0 t , 0) and p A = s/2 (0, 0 t , 1), where s is the squared center of mass energy of the p+A system. The energy (or longitudinal momenta) fractions x 1 and x 2 of the incoming gluons from the projectile and the target, respectively, can be expressed in terms of the rapidities (y 1 , y 2 ) and transverse momenta (p 1t , p 2t ) of the produced particles as where m denotes the quark mass. By imposing production in the forward direction, we effectively select these fractions to be x 1 ∼ 1 and x 2 1. Therefore, the large-x gluons of the dilute projectile are described in terms of the usual gluon distribution of collinear factorization g(x 1 , µ 2 ), and the pA → QQX cross section is obtained from the gA → QQX cross section as: where p = (p + , p t ) denotes the momentum of the incoming gluon. By contrast, due to the large gluon density of the small-x 2 gluons, the gA → QQX cross section does not generally factorize further (it does if non-linear effects can be neglected): dσ(gA → QQX) = dσ(gg → QQX) ⊗ g A . This is due to the fact that in the saturation regime, the gluons in the nuclear wave function interact with the projectile in a coherent manner. Such density effects can be taken into account using the CGC description of the dense small-x 2 gluon content of the nucleus in terms of strong classical fields. Then, the gA → QQX cross section involves averages over color field configurations which may be written as where W x2 [A − ] represents the probability of a given field configuration (we use a gauge in which A − is the only non zero component of the field). Let us now detail what these CGC averages exactly look like, and how an effective factorization with several TMDs emerges in the appropriate limit [42]. . Amplitude for quark-antiquark production in the CGC formalism. The pair can be radiated before (left) or after (right) the interaction with the target. The two terms come with a relative minus sign.

A. Starting CGC formulation
Our starting point is the CGC formalism for quark-antiquark pair production in dilute-dense collisions. The amplitude for quark-antiquark pair production is schematically presented in Fig. 1. In the CGC formalism, the scattering of the partons from the dilute projectile with the dense target is described by Wilson lines that resum multi-gluon exchanges; fundamental Wilson lines for quarks and adjoint Wilson lines for gluons. As a result, the cross section involves multipoint correlators of Wilson lines. In particular, the square of the amplitude from Fig. 1 contains four terms: a correlator of four Wilson lines, S (4) , corresponding to interactions happening after the creation of the qq pair, both in the amplitude and the complex conjugate, then a correlator of two Wilson lines, S (2) representing the case when interactions with the target take place before the gluon splits in both the amplitude and the complex conjugate one, and two correlators of three Wilson lines, S (3) , for the cross terms. Introducing the cross section reads [12]: where denote the transverse positions of the final-state quark in the amplitude and the conjugate amplitude, respectively, and denote the transverse positions of the final-state antiquark in the amplitude and the conjugate amplitude, respectively. The difference u −u is conjugate to the hard momentum P t , and v −v is conjugate to the total transverse momentum of the pair k t . The S (i) Wilson line correlators are given by: where with t a and T a denoting the generators of the fundamental and adjoint representation of SU (N c ), respectively. The functions ϕ λ αβ are the g → QQ splitting wave functions, and their overlap is given by: with m denoting the mass of the quark and with B. Extracting the leading power In order to investigate the TMD regime, we shall extract the leading power in 1/P 2 t . This corresponds the quark and the antiquark being emitted nearly back-to-back in the transverse plane, as the total transverse momentum of the pair |k t | is required to be much smaller than the individual transverse momenta. Importantly, even though Q 2 s is also required to be much smaller than P 2 t , saturation effects still play a role, when k 2 t ∼ Q 2 s . In the case of massless quarks, the calculation was performed in [12] in the large-N c limit and in [22] keeping N c finite.
In the |k t |, Q s |P t | limit, the integrals in (6) where Then, the combination inside the brackets . in Eq. (6) can be rewritten: This expression vanishes if either u or u is set to zero. Therefore, the first non-zero term in its expansion is the one that contains both one power of u and one power of u : So far the derivation has been identical to that of the massless case in [22], the difference resides in the wave function overlap (13), which, after multiplication by u i u j and Fourier transformation, yields: In the massless case, the transverse indices of the various structures in (20) were projected onto δ ij only, and unpolarized gluon TMDs were emerging. Now the presence of the mass is responsible for the appearance of new objects: the so-called linearly-polarized gluon TMDs.

C. Unpolarized and linearly-polarized gluon TMDs
The last integrations which remain to be done correspond to definitions of various gluon TMDs: Both parts of the projection are gluon TMDs, as we will shortly demonstrate. Interestingly, the traceless parts -H -are the TMDs that correspond to the linearly polarized gluons inside the unpolarized nucleus [2,11,14,40]. Gluon polarization hence does play a role in forward heavy-quark production in dilute-dense collisions, even when those collisions involve unpolarized beams. The connection between the generic operator definitions of the gluon TMDs and the definitions given here, valid in the small-x limit, was detailed in [22] (strictly speaking after projecting onto δ ij , but the derivation is identical otherwise), and a short summary can be found in appendix B.
In the leading-logarithmic approximation, the evolution of the CGC wave function In turn, CGC averages in general, and the 6 gluon TMDs introduced above in particular, also evolve towards small values of x 2 according to that non-linear equation. In addition, the scale dependence of those gluon TMDs (not made explicit here) can also be taken into account, although we leave for future work: at small-x this boils down to implementing Sudakov factors into our formalism [43,44].
Introducing the angle φ between P t and k t , one can write and put all the pieces together to finally obtain: Therefore, the leading power of the CGC expression can be interpreted as a TMD factorization formula, with the small-x 2 gluon carrying a transverse momentum equal to k t , and with several gluon TMDs needed to consistently describe the dense gluon content of the nucleus. This is illustrated by Figure 2. As is clear from the above formula, the information on the gluon polarization, encoded in H , couples to the mass m of the heavy quarks, and exhibits an angular dependence cos (2φ), where φ is the angle between the transverse momentum of one of the jets, and the transverse-momentum imbalance of the two jets.  Figure 2. One of the leading order diagrams for inclusive heavy-quark pair production in p+A collisions.

D. Final formula
It is worth noting that F gg , and H (2) gg are not independent, but instead are related to each other through the dipole distribution in the adjoint representation: F ADP (x 2 , k t ) (different from the fundamental dipole gluon TMD F (1) qg = F DP ), defined as: Indeed, we have: as well as: Therefore, the cross section may finally be written as: This is our final formula. It is no more complicated than the one derived in [15] (3 Fs and 2 Hs) (which we show how to recover below), but it is more general (we did not assume the MV model and we kept the complete quark-mass dependence). Moreover, we clearly show that it is the adjoint-dipole gluon distribution which is involved. This TMD is different from the more familiar fundamental-dipole gluon TMD, but it features the same property that, in the small-x limit, its unpolarized and linearly-polarized versions are identical.
In section IV, we shall evaluate numerically all the three unpolarized gluon TMDs as well as their linearly-polarized partners using a lattice calculation to solve the JIMWLK equation We will also evaluate F ADP directly and check that F

A. McLerran-Venugopalan model
The McLerran-Venugopalan (MV) model [45][46][47] (see also [48]) is a classical model for the gluon distribution in a large nucleus. It assumes a Gaussian distribution of color charges, which act as static sources, generating the soft gluons through the Yang-Mills equations. The two-point function of the gluon field A − a is given by: with where λ A (x + ) is the density of color charge squared of the valence quarks, per unit volume and per color. Its precise dependence on x + is not important, since all final results only depend on the integrated density µ A , given by: Evaluating the dipoles, defined in Eqs. (11) and (18), within the MV model yields: and S where we introduced the dimensionless quantity: After regulating the infrared, the integral above can be evaluated to logarithmic accuracy, giving: where we have defined: The above transverse momentum scales are the saturation scales experienced by a quark or a gluon, respectively. These definitions allow to write and S (2) gg (r) = e − r 2 Note that the MV model is purely classical, valid for a large nucleus and for values of x 2 of the order of ∼ 10 −2 . Hence, at this point, there is no evolution, which is why the subscript x 2 in the target averages is omitted.

B. Expressions for the gluon TMDs
Let us start with the Weizsäcker-Williams gluon TMD F gg and its partner H gg . As is clear from their definition in Eqs. (33) and (36), their main ingredient is a double derivative of the quadrupole correlator: The expression for the quadrupole in the MV model was obtained in [12], and reads: with: and Plugging expression Eq. (47) into (46), one obtains: To proceed, we can evaluate the derivative of Γ(v − v ) further: which, depending on the projection of the Lorentz indices, gives: where cos α = k t · (v − v ) / (|k t | · |v − v |). Using the integral representation of the Bessel function of the first kind: and the definition of the saturation scale in the MV model, Eq. (44), one then finally obtains the following expressions for the Weizsäcker-Williams gluon distribution and its polarized partner (in accordance with the literature [10,11,14]): where S ⊥ denotes the transverse area of the nucleus: The calculation of the four other gluon TMDs, F gg , F gg , H gg and H (2) gg , is analogous to the one above. Indeed, once again the main ingredient of the gluon TMDs is a double derivative of a correlator of Wilson lines. This time, it is the correlator of the product of two dipoles, which was calculated in the MV model in [49]: and with F the same as in Eq. (48). The gluon TMDs F (1) gg and H gg , see Eqs. (31), (34), are built from the following structure: which, with the help of Eq. (58), becomes in the MV model: Likewise, F gg and H (2) gg are built from the structure: and one obtains: From these expressions, one can write which allows to further obtain F showing that these exact relations are not spoiled by the MV model assumptions. It is also possible to obtain more explicit expressions, using the following intermediate results: and from which, in combination with Eqs. (52) and (53), one finds: (70)

C. Comparison with the literature
In order to recover the results found in [15] from our cross section (26), one needs to introduce the following auxiliary TMDs xG qq (x 2 , k t ) and xH qq (x 2 , k t ), defined as: From (67)-(70), it is straightforward to derive and to write the unpolarized and linearly polarized part of the cross section (26) in the following way (to leading order in m 2 /P 2 t ): These are the expressions derived in [15]. They are valid in the MV model (without the need to neglect the ln(r) dependence of the saturation scale) only, as they make use of (72), but they are not simpler than the general form we obtained in Section II (the same number of TMDs are involved).

A. Numerical implementation
The JIMWLK evolution equation in rapidity, y = ln (x 0 /x), can be solved at fixed coupling in the small-x regime on a two-dimensional lattice with a Langevin diffusion process of SU (3) matrix variables [30]. The matrix degrees of freedom represent partonic Wilson lines along a light-cone direction and the lattice discretizes transverse space. We use the numerical code developed for the calculation of unpolarized gluon TMDs in [22]. This code is based on algorithms described by Rummukainen and Weigert [50] and Lappi [51,52]. We choose to generate the initial SU (3) configurations at y = 0 in the McLerran-Venugopalan (MV) model wherein analytical calculations of the gluon TMDs have been performed in the previous section.
The gluon distributions listed in Eqs. (31)-(36) are two-point functions defined as products of various traces, which must be evaluated component-wise with respect to the spatial and color indices, in order to express them as scalar convolution products which can be calculated efficiently using a discrete fast Fourier transform algorithm. The continuum derivative ∂ α can be replaced either by a discrete forward or a central difference operator ∇ α . The most convenient expressions for a numerical implementation on a square lattice of size L × L are the following formulas for the unpolarized gluon distributions: and for the linearly polarized distributions: where k α t is the momentum in the lattice Brillouin zone and k α t is either the forward lattice momentum k α t = 2 sin k α t 2 or central lattice momentum k α t = sin k α t , in accordance with the definition of the difference operator ∇ α .
All these lattice gluon distributions have the correct continuum limit and are thus expressed in terms of 2N 4 c complex two-dimensional discrete Fourier transforms.

B. Adjoint sum rules
With the above definitions, the continuum sum rules (28) and (29) remain true on the lattice. The dipole correlator in the adjoint representation is defined as Taking the second-order lattice cross-derivative with respect to v α and v α of the R.H.S. of (76) we recognize in the various terms, after some interchange (ij) ↔ (kl) of dummy color indices and/or dummy spatial indices v ↔ v , the contributions to F (1) (x, k t ) and F (2) (x, k t ) in (74). Hence we have a lattice sum rule which holds in fact configuration by configuration on a finite periodic lattice where the adjoint correlator S This sum rule remains valid whether one defines the lattice derivative as a forward (k The lattice adjoint sum rules hold true configuration by configuration only if one uses the same definition of the lattice derivative in both sides of (78). Since it requires only one additional discrete Fourier transform from a single Langevin simulation, using a different definition provides a very economical way to measure the growth of lattice (and rapidity) discretization errors during the JIMWLK evolution.

C. JIMWLK evolution
All numerical measurements of the gluon distributions studied in this work have been performed on the lattice size L = 1024. A statistical sample of 50 independent trajectories in rapidity has been generated with initial SU (3) configurations distributed so that the correlation length of the dipole correlator in the fundamental representation be in lattice units: This correlation length is defined following the standard Gaussian-like convention: It is indeed unreliable to study lattice momenta beyond π 4 in the Brillouin zone of a periodic lattice and a safe upper bound for the initial correlation length is R s L 16 (the correlation length then decreases with evolution). We treat the issue of those lattice artifacts in momentum space which are due to the breaking of O(2) rotational invariance by a square lattice exactly as explained in [22]. In all subsequent figures, we choose to keep only the (integer) momenta k = (k 1 , k 2 ) with All averages or data points displayed in the plots have error bars which have been determined from the JIMWLK evolution of the random sample with a Langevin step δs = α s π 2 δy = 10 −4 . To interpret the results it is convenient to fix physical units. If we assume a starting x of x 0 = 10 −2 , with associated saturation scale Q 2 s (x 0 ) = 0.2 GeV 2 , we can restore the lattice spacing a from (80) and find: Choosing α s = 0.15, a value of α s /π 2 y = 0.1 would correspond to x 1.4 10 −5 . Fig. 3 displays the initial gluon TMDs calculated on the lattice in the MV model at y = 0. As already observed in [22] for the unpolarized gluon TMDs F (1) gg and F (3) gg , the linearly polarized gluon TMDs H (1) gg and H (3) gg , as well as the adjoint dipole correlator F ADP , have also the expected universal 1/k 2 t behavior at large k t [21]. The high-k t behavior during the JIMWLK evolution is best exhibited by looking at the top plot of Fig. 4, which shows the gluon TMDs for α s y/π 2 = 0.1, after enough evolution to have reached the geometric scaling regime, but not too much so that the high-k t tails of the gluon distributions stays within the accessible momentum range on the lattice. Our results for F (1) gg and F (3) gg match indeed with [22]. Furthermore, we reproduce, at least qualitatively, the numerical results for F (3) gg and H (3) gg in [53]. The important observation to be made is that the data confirm the observations made earlier in [22], that in the limit of large k t , the high-energy or k t -factorization regime is recovered, in which all gluon TMDs converge to a common unintegrated PDF. The only exceptions are the distributions F We also show in the bottom of Fig. 4, the gluon TMDs after further evolution at α s y/π 2 = 0.2, where the high k t has disappeared from the accessible momentum range of our analysis. On the other hand, this allows us to probe the saturation regime at low k t , where the various gluon TMDs are very different from each other, and where the process dependence of TMDs is most relevant and cannot be ignored. Furthermore, when evolving towards smaller values of x, the gluon TMDs shift towards larger values of k t . This is to be expected from the fact that the distributions follow the saturation scale, which grows when x decreases. The different values of α s /π 2 y, for which we plot the numerical results, are listed in Table I, along with the corresponding value of x, as well as the approximate value of the saturation scale Q sg directly measured from the maximum of the adjoint dipole distribution (this method of extraction implies small differences with the values of N c /C F Q s = 1.5 Q s obtained with the definition (81)).
Regarding the information on the gluon polarization, H gg and H (2) gg are small to begin with, but the magnitude of H (3) gg is comparable to that of F gg is equal to F (2) gg and H (1) gg is equal to H (2) gg ; this is a consequence of the adjoint sum rules (78) in combination with the fact that the adjoint dipole correlator F ADP vanishes as k 2 t for k t = 0. The linearly-polarized gluons get suppressed after the first steps in the evolution, however, they are not completely washed out. Indeed, the distributions of linearly-polarized gluons all remain non-zero for momenta of the order of the saturation scale.

V. CONCLUSIONS
In this paper, we used the CGC framework to compute the cross section for the forward production of a heavy quark-antiquark pair in proton-nucleus collisions. In the correlation limit, in which the outgoing quarks are almost back-to-back in the transverse plane, our result could be cast into a TMD factorization formula, involving six different gluon TMDs. Three of these TMDs are unpolarized, and also appear in the cross section for forward dijet production. They are each accompanied by a partner which couples via the quark mass, and which corresponds to the linearlypolarized gluons inside the unpolarized nucleus. We have obtained analytical expressions for each of the TMDs in the MV model. Furthermore, the gluon TMDs were numerically evolved in rapidity using the nonlinear JIMWLK evolution equation.
Our results indicate that the various distributions of linearly-polarized gluons always remain non-zero for values of the gluon transverse momentum of the order of the saturation scale. This observation provides us with a novel way to test parton saturation at the LHC, and to extract the poorly known linearly-polarized gluon distributions. We note that the LHCb detector would be particularly well suited to perform such a measurement of heavy mesons in the forward region, although a detailed feasibility study remains to be done. Indeed, the linearly-polarized gluons impact the production cross section via an angular modulation whose magnitude we plan to better quantify.
In appendix A, we give an outline of the derivation for a similar but simpler process, γ * A → QQX, in which only two of the gluon TMDs appear: the Weizsäcker-Willams distribution F (3) gg (x, k t ) and its polarized partner H (3) gg (x, k t ). Interestingly, the dependence on H gg (x, k t ) via the azimuthal angle between P t and k t , couples not only to the quark mass but also to the virtuality of the photon. This provides alternatives to the pA → QQX process -namely dijets at an Electron-Ion Collider [54], or heavy QQ pair production in ultra-peripheral collisions of heavy ions -using processes whose theoretical formulation involves less gluon TMDs, but which may be experimentally more challenging or more distant in the future.
Finally, let us stress again that the focus of this work was on the implementation of the small-x JIMWLK evolution, and that we did not discuss the scale evolution. That aspect was recently studied in the simpler context of the pA → γ * qX process [55], and it was found that the scale evolution leads to a Sudakov suppression of the angular modulation induced by the linear polarization of gluons. However, in that process only the fundamental-dipole gluon TMD is involved, which is a peculiar TMD since, as we already pointed out, the unpolarized and linearly-polarized distributions are identical at small-x. We leave it for future work to estimate the effect of the scale evolution on the TMDs displayed in Fig. 4, involved in the pA → QQX process. We also note that, at next-to-leading order, additional gluon TMDs appear [56] in the pA → γ * qX process, related to those discussed in this work.
The cross sections (A5) and (A6) were first obtained in [11]. Recently, the next-to-leading power was also obtained in the massless limit [57].

Appendix B: Transverse momentum dependent gluon distributions
In this short paragraph, we demonstrate the equivalence of the small-x Weizsäcker-Williams gluon distribution (for a left-moving hadron) as defined in Eq. (33), and its standard operator definition [4,12]: with Setting exp ix 2 p − A ξ + ≈ 1 since x 2 is small, as well as making use of translational invariance and the fact that the light-cone Fock states are related to the hadronic states |A as follows (see [12,22]): with the normalization A|A = (2π) 3 2p − A δ (3) ( 0), we obtain: Then, using the rules for the decomposition of Wilson lines, such as: U (+∞, w + , w) = U (+∞, −∞, w)U (−∞, w + , w) , the average in Eq. (B5) becomes and one recovers expression (33). The proof for the other gluon TMDs is analogous.