Probing gluon saturation with next-to-leading order photon production at central rapidities in proton-nucleus collisions

We compute the cross section for photons emitted from sea quarks in proton-nucleus collisions at collider energies. The computation is performed within the dilute-dense kinematics of the Color Glass Condensate (CGC) effective field theory. Albeit the result obtained is formally at next-to-leading order in the CGC power counting, it provides the dominant contribution for central rapidities. We observe that the inclusive photon cross section is proportional to all-twist Wilson line correlators in the nucleus. These correlators also appear in quark-pair production; unlike the latter, photon production is insensitive to hadronization uncertainties and therefore more sensitive to multi-parton correlations in the gluon saturation regime of QCD. We demonstrate that $k_\perp$ and collinear factorized expressions for inclusive photon production are obtained as leading twist approximations to our result. In particular, the collinearly factorized expression is directly sensitive to the nuclear gluon distribution at small $x$. Other results of interest include the realization of the Low-Burnett-Kroll soft photon theorem in the CGC framework and a comparative study of how the photon amplitude is obtained in Lorenz and light-cone gauges.


Introduction
Photons produced in proton-nucleus (p+A) collisions are powerful probes of the fundamental many-body structure of strongly interacting matter. Sufficiently energetic photons, in particular, are free of the uncertainties that attend the fragmentation of partons into hadrons. Further, the information on microscopic dynamics that is probed by photon final states complements in fundamental ways the information accessible in deeply inelastic scattering (DIS) experiments off hadrons and nuclei. A global analysis of photon production in DIS and p+A collisions, analogous to those performed for currently for parton distribution functions has therefore the potential to reveal universal features of parton dynamics that are distinct from those that are particular to the scattering process.
The photons produced in deuteron-nucleus collisions at the ultrarelativistic energies of the Relativistic Heavy Ion Collider (RHIC) and p+A collisions at the Large Hadron Collider (LHC) are uniquely sensitive to strongly correlated states of gluons in which the gluons have the maximal occupancy allowed by QCD. The dynamics of these saturated gluon states is described by the Color Glass Condensate (CGC) [1,2,3], a classical effective theory of QCD in the high energy asymptotics (of high center of mass energies and small values of parton momentum fractions x) that may be applicable for significant kinematic windows at RHIC and LHC.
A characteristic feature of the CGC is that the dynamics of saturated gluons is governed by an emergent "saturation scale" Q S (x), which grows with increasing energy, or equivalently, decreasing x. Modes in hadron/nuclear wavefunctions with k ⊥ < Q S are maximally occupied with an occupancy that is parametrically 1/α S , where α S is the QCD fine structure constant. In contrast, the modes k ⊥ Q S have small occupancies and interact dynamically as the partons of perturbative QCD (pQCD). If Q S is large relative f q γ A A Figure 1: Leading order process for prompt photon production in proton-nucleus collisions. The diagram describes the bremsstrahlung of a photon from a valence quark after multiple scattering off the classical gluon field in the nucleus.
to the intrinsic non-perturbative QCD scale, the strongly correlated dynamics of gluons can be computed using weak coupling techniques. In nuclei, the saturation scale (Q A S ) 2 ∼ A 1/3 . The high energy dynamics of nuclei are therefore well suited to test the CGC description of high energy QCD. Photon production, as noted, is a particularly sensitive probe because it is independent of the details of how partons fragment into hadrons.
Photon production in p+A collisions was previously computed to leading order in the CGC framework [4]. The power counting for proton-nucleus collisions corresponds to a "dilute-dense" limit, where contributions to lowest order in Q p S /k 1⊥ 1 (where Q p S is the saturation scale in the projectile proton) are preserved along with all order terms in Q A S /k 1⊥ . Here k 1⊥ , k 2⊥ respectively correspond to the momentum exchange from the proton and the nucleus to the final state of interest. For photon production, the leading term in this dilute-dense power counting corresponds to order O(α), where α is the QED fine structure constant. This leading order (LO) contribution is illustrated in Fig. 1. The quark line in this figure corresponds to a valence quark in the wavefunction of the projectile proton. In the high occupancy regime of k ⊥ ≤ Q A S , there is no α S dependence at LO because the α S factor in the cross section arising from the coupling of a gluon to the valence quark is compensated by the 1/α S occupancy of these gluons in the target. The O(α) dependence must be understood as being accompanied by the valence quark distribution f q in the proton.
At next-to-leading order (NLO) O(αα S ), there are a number of contributions which can be classified into the three classes shown in Fig. 2. The leftmost diagram (class I) corresponds to a gluon emitted from the incoming valence quark. For inclusive photon production, one has to integrate over the phase space corresponding to the emitted gluon. Another NLO contribution to inclusive photon production in this diagrammatic class (not shown) arises from the interference between the LO contribution in the amplitude and a contribution in the complex conjugate amplitude corresponding to the virtual emission and absorption of a gluon by the valence quark. Both of these diagrams give logarithms that are sensitive to the transverse momentum and x of the gluon. These diagrams contribute to the DGLAP renormalization group (RG) evolution of the valence quark distribution 1 . The NLO contributions of class I are therefore actually of O(α) if the bare valence quark distribution in the proton is replaced by the RG evolved quark distribution.
The class II NLO diagram in Fig. 2 was computed recently [5]. In this case, since the quark-antiquark pair are emitted from the gluon prior to their subsequent annihilation into the photon, this diagram is of order O(αα S ) with the cross section for this contribution accompanied by a factor f g corresponding to the gluon distribution in the proton. While this is an NLO diagram, it can in principle provide a much larger contribution to photon production. This is because the gluon distribution grows rapidly while the valence quark distribution decreases at small x, giving f g f q . Thus for photon production in the small x kinematics of the projectile proton, such NLO contributions can dominate significantly. The kinematics of the class II f q γ g Figure 2: Next-to-leading order processes for prompt photon production in proton-nucleus collisions. The class I diagram describes the bremsstrahlung of a photon from a valence quark. The next two diagrams correspond to a gluon from the proton splitting into a qq that annihilates into a photon final state (class II), or emits a photon either before or after rescattering off the nucleus (class III). As described in the text, the class I diagram, upon evolution, is the same order as the diagram in Fig. 1.
diagram is relevant for inclusive photon production at central rapidities. Because of pair annihilation, the transverse momentum of the photon for this diagram is strongly constrained to be dominated by momenta around k ⊥ ∼ Q A S . Thus this contribution is in principle very sensitive to the saturation scale of the nucleus. However as also noted in Ref. [5], there is a further class of NLO processes, class III in Fig. 2, that contribute significantly to photon production in p+A collisions. Firstly, like the class II process, this NLO contribution comes accompanied by a factor f g which, as noted, will overwhelm the LO contribution at small x. Secondly there are some features of the class III computation that are qualitatively different from those of class II. Unlike the latter, photon production, while sensitive to Q A S , is dominated by soft momenta with k ⊥ < Q A S . Similarly, since photon production is not as kinematically constrained for class III diagrams relative to class II diagrams, it will also dominate at large k ⊥ > Q A S . In particular, class III diagrams will match the contribution from leading twist pQCD at high k ⊥ , while the class II contribution is proportional to a higher twist four point correlator even at large k ⊥ . The sum of class II and class III diagrams constitute the relevant NLO contribution to inclusive photon production in the CGC framework.
In this paper, we will compute the class III NLO diagrams for photon production in p+A collisions in the CGC framework 2 . We will perform the computation first in Lorenz gauge ∂ µ A µ = 0 gauge, and subsequently in light-cone gauge A + = 0. In addition to being an independent non-trivial check of our results, the intermediate steps are interesting and realized differently in the two gauges.
The paper is organized as follows. In Section 2, after a preliminary discussion of dilute-dense collisions in the CGC framework, we will outline the derivation of the amplitude for the inclusive production of a qqγ final state in Lorenz gauge. Our work closely follows the previous derivation in this gauge of the amplitude for gluon production [7] and quark-antiquark pair production [8] in proton-nucleus collisions. As for the case of the pair production amplitude considered previously, we show that contributions from so-called "singular" terms, wherein the photon is produced from within the target, are exactly canceled by terms that one can identify as gauge artifacts in regular terms. The latter are contributions where the quark-antiquark pair (and the photon) are produced either before or after the scattering of gluons from the target off the projectile.
In Section 3, we compute the cross section for inclusive photon production. For readers interested in the central result of this work, the key expression is given in Eq. (62). As in the case of quark-antiquark production [8], the cross section factorizes into the product of the unintegrated gluon distribution in the projectile times the sum of terms corresponding to the gluon distribution in the target, and quark-antiquarkgluon and quark-antiquark-quark-antiquark light-like Wilson line correlators. These correlators contain nontrivial information on many-body gluon and sea quark correlations that are of all twist order in conventional pQCD language. In Sec. 4 we demonstrate that for Q A S k 2⊥ , the cross section can be expressed as a k ⊥ factorized product of unintegrated gluon distributions in the projectile and target [9,10]. We show explicitly that the corresponding leading twist k ⊥ factorized amplitude agrees exactly with the expression derived recently by Motyka, Sadzikowski and Stebel [10] in the context of Z 0 boson hadroproduction. Taking the limit k 1⊥ , k 2⊥ → 0, one recovers the formal structure of the leading collinear factorization contribution to inclusive photon production from gluon-gluon scattering [11]. This demonstrates the unique sensitivity of inclusive photon production to the nuclear gluon distribution function at small x. Unfortunately, the detailed analytical comparison of the two expressions is cumbersome; a numerical comparison is left for future work. In Sec. 5, we summarize our results and provide a detailed outline of the necessary ingredients for the numerical computation of inclusive photon production in proton-nucleus collisions. We will leave this numerical computation for a subsequent publication.
In Appendix A we will compute the amplitude for inclusive photon production in light-cone gauge A + = 0. This computation follows a previous computation of gluon production in this gauge [12]. Unlike gluon production, where the amplitudes in Lorenz gauge and light-cone gauge differ (albeit the cross sections of course agree), in this case, the amplitudes in the two gauges agree exactly. Some properties of the amplitude are outlined in Appendix B. In addition to demonstrating that this amplitude satisfies a Ward identitiy, we show in particular that the well known Low-Burnett-Kroll theorem [13,14,15] is satisfied. This corresponds, in the soft photon limit, to the factorization of the amplitude into the amplitude for quark-antiquark production times a contribution determined by the Lorentz structure of the photon.

Amplitude for inclusive qqγ production in the Lorenz gauge
In the CGC effective theory, the gluon dynamics with high occupancy at small x are described by the classical Yang-Mills equations. For collisions of the proton moving in the positive z direction and the nucleus moving in the negative z direction at the speed of light, the Yang-Mills equations are where ρ p and ρ A correspond to the static and stochastic color charge density matrices at large x in the proton and the nucleus, respectively. This framework is feasible at a given impact parameter for three different regimes of color charge densities depending on the energy of the collision and the kinematic regime of interest. These correspond to the dilute-dilute (ρ p /k 2 1), the dilute-dense (ρ p /k 2 1⊥ 1, ρ A /k 2 2⊥ ∼ 1), and the dense-dense (ρ p /k 2 1⊥ ∼ 1, ρ A /k 2 2⊥ ∼ 1) regimes. To compute inclusive cross sections, the modulus squared of the computed amplitude must be averaged over the color sources, thereby replacing ρ p , ρ A → Q p S , Q A S in the power counting. Subsequent to this averaging, the dilute-dilute limit in the CGC effective theory can be matched to the leading twist frameworks of pQCD at small x. Further, the dilutedense limit corresponds to the dominance of leading twist contributions on the proton side with all twist contributions on the nuclear side, and the dense-dense limit includes all twist contributions from both the proton and the nucleus to the scattering process of interest.
At the collider energies currently accessible in p+A collisions, the appropriate dynamics is that of the dilute-dense limit. The nucleus is dense because there is an enhancement of the order of A 1/3 in the number of color sources even at large x. These sources can all radiate soft gluons that further function as color sources for even softer gluons. In the proton, on the other hand, there are O(1) color sources, and it is only at very large rapidities relative to the proton fragmentation region that the density of color sources becomes ρ p /k 2 1⊥ ∼ 1. Alternately, the dense-dense regime can be attained in very rare high multiplicity events where the proton fluctuates into sources of color charge already at large x.
In this work, we will be interested in photon production in p+A collisions at collider energies for the small x kinematics, and in event classes where the dilute-dense asymptotics is appropriate and gives the dominant contribution. Analytic computations are feasible then, and explicit expressions can be derived. In contrast, the dense-dense limit is not analytically tractable even for inclusive gluon production since there is no small parameter to expand in. Results can be obtained only through numerical simulations of the Yang-Mills equations [16,17,18].

Structure of the gluon field and setup of the amplitude computation
In what follows, we work entirely in the Lorenz covariant gauge ∂ µ A µ = 0. The outline of the derivation in light-cone gauge is given in Appendix A. For a dilute-dense system with ρ p /k 2 1⊥ 1 and ρ A /k 2 2⊥ ∼ 1, the analytical solution of Eq. (1) in the Lorenz gauge can be expressed in momentum space as [7] A µ (q) = A µ p (q) + ig where A µ (q) = A a µ (q)T a and T a are the generators of SU(N c ) in the adjoint representation. We will also use a shorthand notation, k ⊥ ≡ d 2 k ⊥ (2π) 2 and x ⊥ ≡ d 2 x ⊥ , hereafter. In detail, the ingredients going into Eq. (2) are as follows. The first term, A µ p (q), represents the gluon field of the proton alone The vectors C µ U and C µ V are abbreviated forms for the momentum dependence of the integral, with k 1⊥ being the momentum exchanged from the proton, while k 2⊥ ≡ q ⊥ − k 1⊥ exchanged from the nucleus. The explicit forms of C µ U and C µ V are These two functions are related to the well-known Lipatov effective vertex [19,20,21] via the simple relation, This effective vertex is a gauge covariant expression that efficiently combines the various contributions to glue-glue scattering in the Regge-Gribov limit of QCD.
In Eq. (2), the Wilson lines U and V account for the modification of the gluon field of the proton due to multiple gluon scatterings in the nucleus and can be expressed, for an arbitrary light-like path as where A − A (z + , x ⊥ ) represents the gluon field alone. This field can be expressed similarly to A µ p above as In the above expressions for the Wilson lines P + denotes path-ordering in the x + direction. While U is the standard Wilson line in the adjoint representation of SU(N c ), V is an unusual form of the Wilson line that turns out to be a gauge artifact of the gluon production amplitude in Lorenz gauge. This can be deduced from the fact that the V dependent terms in the amplitude do not appear in the cross sections nor in calculations in other gauges [12,22]. While one expects therefore that V should drop out at the level of the cross section for photon production, we will show that they will not appear in our final expressions for the amplitude. The same observation was made previously for the amplitude for qq production in Lorenz gauge. We will henceforth use the following shorthand notation for the complete Wilson lines, Before we proceed, it is instructive to discuss the underlying structure of Eq. (2) in coordinate space. This can be decomposed by splitting Eq. (2) as A µ = A µ R + A µ S where A µ R is called a regular term and A µ S k 1 : 4-momentum exchanged from the proton k 2 : 4-momentum exchanged from the nucleus k γ : photon 4-momentum k: 4-momentum fromŨ q: quark 4-momentum p: antiquark 4-momentum P : 4-momentum of the final state P = k 1 + k 2 = p + q + k γ (4-momentum fromŨ † : a singular term. The former corresponds to the emission of the gluon from the proton before interacting with the highly Lorentz contracted shock wave of gluons that comprises the nuclear target, while the latter corresponds to gluon production from within the shock wave. The latter term is proportional to the Lorentz contracted width of the nucleus ∝ δ(x + ). Because 2q − = q 2 /q + + q 2 ⊥ /q + , we can split C − V (q), the only term with q − dependence, into a regular part C µ V,reg (q) and a "singular" part We observe that the singular field appears from the second term in Eq. (9) because the q 2 term in Eq. (9) cancels the 1/q 2 pole in Eq. (2). Clearly, the q − integration leads to δ(x + ), but we need to be cautious about the infinitesimal longitudinal extension in V (x ⊥ ). After careful treatment, its Fourier transformed representation in coordinate space can be expressed as In addition to decomposing the gauge field, we will need one more ingredient for our computation of the class III amplitude for photon production. This is the effective vertex corresponding to self-energy insertions in the time-ordered quark propagator arising from multiple insertions of the nuclear gluon background field: it is given by [8] T where p is the quark momentum and k is the momentum transfer from the multiple gluon "kicks" to the quark. The Wilson line in the fundamental representationŨ is defined as in Eq. (5) with the SU(N c ) adjoint generators T a replaced by the generators t a of the fundamental representation. Explicitly, it is written as [23] U ( The effective vertex appears in the quark propagator andŨ andŨ † sum over all the gluon insertions to the quark and the antiquark respectively. The momentum transfer from the nucleus k should be integrated over. Since the effective vertex does not change the order of the diagrams parametrically by any power of g, processes where the emitted photon is sandwiched between two such effective vertices are possible. We will see in the following that this is not the case; the kinematics of the process will constrain the number and configuration of these effective vertices.
From the above discussion, we can deduce that for the amplitude for photon production in Fig. 2 there are fourteen non-vanishing contributions. We will separate them in two groups of diagrams: 1. the regular terms, in which the qq pair is created from the regular field A µ R . In the following computation, we shall include only one insertion of A µ R/S on the quark propagator, regular or singular, to stay consistently at first order in the proton source ρ p . As noted, the number of gluon insertions from the nucleus onto the quark and antiquark lines from the nucleus can, however, in our dilute-dense power counting, be as many as the kinematics permit. The amplitude can be decomposed into the external polarization vector of the photon and an amplitude vector, with the results in the following subsections expressed in terms of the amplitude vector defined as

the singular terms, wherein the pair is spawned from
where q, p, and k γ are the quark, the antiquark, and the photon external three momenta, respectively, and λ is the photon polarization. We summarize in Table 1 the momenta notations that will be used in the following calculations.

Regular contributions to the amplitude
Following the above stated classification, we will proceed to find the regular diagrams which have no fundamental Wilson lines first. For this case, there are two diagrams, shown on Fig. 3, with exactly one insertion of the proton source, in the form of the regular field A R . The two diagrams represent the scattering of a gluon off the target and the resulting creation of a qq pair. The photon is then emitted from the quark or antiquark line. Using standard Feynman rules, the vector amplitude for the diagram (R1) (denoted as where P is the total external 4-momentum and the quark lines are given in this calculation by the vacuum time-ordered fermion propagator The proton field A µ p does not contribute to this diagram. It contains the delta function δ(p − + q − + k − γ ) which cannot be satisfied if the quark, antiquark, and photon are on-shell, as p − , q − , k − γ > 0. Dropping the A µ p term, we are left with the rest of A µ R , which, for the amplitude M µ R1 , gives, A RQ Following the same procedure as the one for (R1), one finds the amplitude contribution (R2) for the photon emitted from the antiquark to be It will be shown in the next subsection that the singular contributions will cancel the terms with the Wilson line V -the final expression has no dependence on V .
Diagrams with one insertion of the effective vertex on the quark propagator can have one Wilson lineŨ orŨ † in the fundamental representation, as shown in Fig. 4, for the quark [diagrams (3)-(5)] and likewise for the antiquark [diagrams (6)- (8)]. In the following steps we will treat them separately for convenience. As in the case of the amplitudes (1) and (2), the amplitude (3) in Fig. 4 will have a regular field insertion. However, now in addition we must insert the effective nuclear vertex (11) for the multiple gluon scatterings. We should integrate over nuclear momentum transfer k 2 to obtain, We then integrate over k + 2 and k − 2 . This integration is trivial for k + 2 since T (q + k γ , k 2 ) contains δ(k + 2 ). Only the proton field part A µ p of the regular field gives a finite contribution. In this part, the k − 2 integration is also trivial because A µ p contains δ(P − − k − 2 ). We shall now demonstrate that the integration over the remaining part of A µ R vanishes by the residue theorem. The singularities in k − 2 of the regular field and the quark propagator S 0 (q + k γ − k 2 ), respectively, are where the first pole comes from the prefactor and the second from C + U in A µ R . Since all the three poles are above the real k − 2 axis, the k − 2 integral vanishes by the residue theorem. This result has the clean physical interpretation that the quark-antiquark pair is first created in the proton and subsequently the quark scatters off the gluons in the nucleus. The calculation of the diagrams (4) and (5) is analogous. We can collect the results for the amplitude vector in a compact form as where the integration variable was changed from k 2 to k 1 = P − k 2 and β ∈ 3, 4, 5. The symbol R µ β (k 1⊥ ) is a shorthand notation for the Dirac structure where we define / p ⊥ ≡ p i γ i and M 2 (p ⊥ ) ≡ p 2 ⊥ + m 2 as the transverse mass. The reductions to transverse Dirac matrices in the numerators comes from using k + 1 = P + , k − 1 = 0 and also the identities (γ + ) 2 = 0, (γ − ) 2 = 0.
The next three contributions to the photon amplitude, represented in diagrams (6)-(8) of Fig. 4, can be found in the same fashion, and as the former, can be shown to have vanishing pole integrations. The amplitudes for these processes are where β ∈ 6, 7, 8 corresponds to the respective diagrams in Fig. 4 and the corresponding Dirac structures To tidy up our notation, we will express the sum of the contributions (3)-(8) in Fig. 4 as where the total Dirac structure is combined as In the first term in (25), we introduced a dummy integration over y ⊥ and k ⊥ . In the second term in (25), we renamed x ⊥ → y ⊥ and further, introduced a dummy integration over x ⊥ and k ⊥ . Let us now consider the case where there are two insertions of the effective vertex on the quark propagator. The contribution corresponding to a photon emission between two insertions of the effective vanishes for the same kinematic reasons as previously -the pole integration yields a null contribution. Thus the only non-zero contributions come from diagrams where one insertion is on the quark line and the other on the anti-quark line. There are four such contributions, which are listed as diagrams (9)-(12) in Fig. 5.
All of these diagrams are computed with the same logic as previously. As an example, we focus on diagram (9). The corresponding amplitude can be written as As for the case with only one effective vertex, and for the same reasons articulated there, only the proton piece A µ p of the regular field A µ R contributes. The integrals over k + , k + 1 and k − 1 can be performed thanks to the δ-functions in the two effective vertices T and in A µ p , respectively. The remaining integration over k − can be evaluated by the method of residues. As in the case of the one effective vertex insertion, the gluon scattering vertices (C µ U and C µ V ) are kinematically forbidden by the pole integration. Performing similar steps for the remaining diagrams (10)- (12), the amplitude of the sum of the diagrams (9)-(12) can be expressed as where The Dirac structures here for β = 9, . . . , 12, corresponding to the diagrams (9)-(12) in Fig. 5, are For clarity of presentation, the following functions in the denominator have been defined as In all the equations presented here, k ⊥ stands for the momentum transferred from the dense nucleus to the quark line. Likewise, k 1⊥ is the momentum transferred from the proton. This is transparent since the proton color sources ρ a p (k 1⊥ ) are dependent on this integration variable alone. From this fact, and from the fact that the initial momentum flow must be inferred from the final state momenta P ⊥ , one can readily notice that P ⊥ − k ⊥ − k 1⊥ is the momentum transfer from the nucleus to the antiquark.

Singular contributions to the amplitude
The terms for which the qq is produced by the singular part of the field A µ S are represented by the diagrams (S1) and (S2) in Fig. 6, corresponding to both the creation of the qq pair, and the subsequent emission of the photon from within the nucleus. The amplitude for this process is non-vanishing in the Lorenz gauge ∂ µ A µ = 0 and needs further explanation. Without regularization, the expressions for amplitudes (S1) and (S2) would look the same as in Eq. (14), but with the regular field exchanged for the singular field. With regularization, the δ(x + ) function has a small width: δ(x + ) → δ ε (x + ). This allows the qq to undergo multiple gluon scatterings.
Diagrams (S1) and (S2) split into four different parts; each one corresponds to insertions on the quark and the antiquark lines, making four distinct combinations. However these insertions have to be treated differently than in the previous terms we considered as they occur inside the regularized region. This can be achieved by changing the Wilson lines in the insertions into incomplete Wilson lines. Summing all the terms that make up (S1), we find, Several points have to be noted about this expression. Firstly, the structure of the amplitude reflects the fact that for these diagrams the quarks do not rescatter again as a consequence of the vanishing duration of the interaction in the ε → 0 limit. In this limit the insertion of the singular field and the rescattering occur in the same transverse plane, which is intuitively what one would expect from a qq pair being created and interacting inside a heavily boosted nucleus. Secondly, the photon emission from inside the nucleus or followed by another scattering would be tantamount to resolving the nuclear gluon shock wave; this too is not kinematically viable. Therefore one can only have the addition of an external photon leg emitted by the outgoing quark or antiquark. Using the identityŨ and the formula [8] ig 2 we obtain, A similar procedure for (S2) gives In these expressions, we have added and substracted the adjoint representation identity matrix to show their similarity to Eqs. (17) and (18). In the following, it will be shown that some of these contributions cancel out.

Assembling the contributions: the complete result for the photon amplitude
The net contribution to the amplitude is given by a summation of the terms in Eqs. (17), (18), (25), (28), (35) and (36). There are several cancellations that occur when we put these terms together. The cancellation of the Wilson line V (x ⊥ ) will be addressed first, as it was anticipated, and it stands as a good check of our computation. Firstly, using the relations of Eqs. (4) and (9) C µ V,reg (P ) is explicitly written as Using the identitiesū it follows that the first term of (37) cancels out within Eqs. (17) and (18). These relations leave us with only the second term in C µ V,reg (P ). It is straightforward to check that both remaining contributions in Eqs. (17) and (18) are the counterparts of the V dependent terms in the singular diagrams, and so Eqs. (35) and (36) demonstrate the anticipated cancellation.
The resulting effective vertex for the U terms in the sum of the diagrams (R1), (S1) and separately (R2), (S2) is which can equivalently be written as where the expression C µ L is the Lipatov effective vertex The terms that survive the above mentioned cancellation give the net amplitude for a gluon to first scatter off the nucleus before emitting a qq pair and can be expressed as with β = 1, 2 and where Combining all the results, the full amplitude vector is with This last expression extends the Dirac structure found in Ref. [8] to photon production. We note that in Eq. (44) we introduced a dummy integration over x ⊥ , y ⊥ and k ⊥ to ensure that all the terms have identical integration variables. The result (44) can be further simplified by making use of the identities, Using Eq. (46), the expression for the amplitude can considerably simplify to This expression for the photon amplitude is a key result of this work. In Appendix A, we will show that this expression for the amplitude in Lorenz gauge is identical to the expression derived in light-cone gauge. Before we conclude this section, we wish to make a few points regarding the final result. Firstly, in Eq. (46), the sum of the four effective vertices is zero as a consequence of momentum conservation. If there is no nuclear and proton momentum transfer, the quark-antiquark dipole cannot be created. The second and third relations in Eq. (46) stand for the vanishing of contributions if there is no momentum transfer either from the projectile or the target. We are thus left only with contributions to the photon amplitude that have i) both the quark and the antiquark interact with the nucleus after being created (and before or after radiating the photon), and those ii) where the gluon from the proton scatters off the nucleus before creating the quark-antiquark pair and thence, the photon.

The inclusive photon cross section
The probability for creating a qq pair with 4-momenta q and p, respectively, and a photon with momentum k γ , for a fixed distribution of sources ρ p in the projectile and ρ A in the target, respectively, is given by Here E p , E q and E kγ denote the relativistic energies of the antiquark, quark and photon, respectively. The sum over polarizations can be taken by noting that The color average of an inclusive quantity O must be taken after taking the modulus squared of the amplitude [7], The weight functionals W p [ρ p ] and W A [ρ A ] are density matrices that obey the JIMWLK evolution equations [24,25,26,27] that describe the renormalization group evolution of distributions of color charges in the wave-functions of the projectile and the target, respectively, from their respective fragmentation regions at large x down to the small x values probed by measurements in high energy collisions.
Thus the impact parameter-dependent triple differential probability corresponding to the probability functional defined above is where d 6 K ⊥ ≡ d 2 p ⊥ d 2 q ⊥ d 2 k γ⊥ and d 3 η K ≡ dη p dη q dη kγ . The angle brackets, · · · , represent the color average in Eq. (50). The triple-differential inclusive cross section can be obtained by integrating the above expression over impact parameter b ⊥ , Using Eqs. (47) and (51), the triple differential probability can be expressed as where we wrote four Dirac traces in a compact form using the following notation, with n, m ∈ {g, qq} and where we used the abbreviations T g ≡ T µ g (k 1⊥ ), T µ g ≡ T µ g (k 1⊥ ) and T µ qq ≡ T µ qq (k 1⊥ , k ⊥ ), T µ qq ≡ T µ qq (k 1⊥ , k ⊥ ). We will now shift the spatial transverse coordinates by the impact parameter to make them relative to the center of the nucleus, Since the Wilson line correlators are approximately translational invariant for a large nucleus, the phases in the integrand yield a factor of e ib ⊥ ·(k 1⊥ −k 1⊥ ) after the above coordinate shift. Thus, to obtain the tripledifferential cross section we can easily integrate over the impact parameter b ⊥ , resulting in a δ-function: (2π) 2 δ (2) (k 1⊥ − k 1⊥ ). In Eq. (53), we also introduce a dummy integration over k 2⊥ , together with a δ-function representing overall momentum conservation (2π) 2 δ (2) (P ⊥ − k 1⊥ − k 2⊥ ). As discussed previously in Refs. [7,8], the color averages in the formula above can be re-expressed in terms of novel unintegrated distribution functions. The correlator of color sources in the projectile proton is defined to be where dϕ p /d 2 r ⊥ is the proton unintegrated gluon distribution per unit area, r ⊥ is the variable that runs over the transverse profile of the proton, C F ≡ (N 2 c − 1)/2N c is the SU(N c ) quadratic Casimir in the fundamental representation. When the momenta in the argument of dϕ p /d 2 r ⊥ are identical as is the case after the integration over b ⊥ , the expression greatly simplifies to with ϕ p (k 1⊥ ) being the unintegrated gluon distribution in the proton. This identification of the correlator of color sources can be straightforwardly generalized to correlators of Wilson lines which likewise can be expressed in terms of nuclear unintegrated distribution functions [8,28,29]. In conventional pQCD language, these unintegrated momentum distributions resum a sub-class of all twist correlations in the nucleus. The correlator of two adjoint Wilson lines can be expressed as Similarly, the three point fundamental-adjoint Wilson line correlator can be expressed as and likewise for its Hermitean conjugate expression in the cross section. Finally, the four point correlator of fundamental Wilson lines can be expressed as Using these relations, the triple-differential cross section can be finally expressed as We remind the reader that P ⊥ = p ⊥ + q ⊥ + k γ⊥ . It is worth noting that the r. h. s. dependence on the external momenta lies only in the overall δ-function and in the Dirac traces. The expression in Eq. (61) would be useful if it were feasible to measure direct photons in coincidence with a quarkonium state, for instance the J/Ψ meson. Alternately, by integrating over either the quark (antiquark), this cross section would provide the rate for direct photons measured in coincidence with open charm (anticharm) states. These measurements are challenging even at LHC energies. On the other hand, inclusive prompt photon differential cross section has already been measured at RHIC in deuteron-gold collisions [30,31] and photon measurements can be anticipated at both RHIC and LHC in p+A collisions in the near future.
The simplest quantity that can be measured is the inclusive prompt-photon single-differential cross section. For this one must integrate over the quark and the antiquark momenta and rapidities, and one obtains, This expression is the main result of this paper. In the previous related work on quark-antiquark production [8,29], the integration over the momentum of the quark or the antiquark enabled one to simplify the result so that it depended only on the φ g,g A and φ qq,g A nuclear distributions. Unfortunately, it appears no such simplification is possible here due to the particular momentum dependence of the τ n,m functions. Equation (62), however, simplifies considerably in the k ⊥ factorization and collinear pQCD limits, which will be the subject of the next section.

k ⊥ factorization and collinear factorization limits
Our master expression, Eq. (62), can be simplified considerably in a high transverse momentum expansion of the nuclear unintegrated distributions. The high momentum expansion corresponds to expanding the Wilson lines,Ũ (x ⊥ ) and U (x ⊥ ) in the fundamental and the adjoint representations, to lowest non-trivial order in terms of ρ A /∇ 2 ⊥ . This is equivalent to the leading twist expansion in pQCD as will become manifest when we consider the collinear factorization limit of our expressions. Physically, this limit corresponds to the dynamics when the density of color sources is large enough to be represented as classical color charges, but only one of the color charges in this classical color distribution is resolved by a high transverse momentum probe.
A further simplification occurs for large nuclei with Gaussian random color sources [32,33,34] that satisfy the relation, where the color charge squared per unit area is µ 2 A = A/(2πR 2 ) ∼ A 1/3 . In general, the correlations amongst color sources in the nuclear wave-function is given by the weight functional W A [ρ A ], which as we noted previously, satisfies the JIMWLK equation. This equation is formally equivalent to the Balitsky-JIMWLK hierarchy [35,36] of n-point Wilson line correlators. The JIMWLK equation has been solved numerically [37,38]. It was found, to a good approximation, that the solution is well represented by a non-local Gaussian [38], with µ 2 such that ϕ A (Y A , k ⊥ ) obeys the Balitsky-Kovchegov (BK) equation [35,39]. In this expression, Y A is the rapidity, relative to the nuclear beam rapidity, of gluons from the target off which the qq pair scatters. The BK equation, for large N c , is the closed form equation for the "dipole" correlator of Wilson lines -the lowest term in the Balitsky-JIMWLK hierarchy. Even though the BK equation has a closed form, it cannot be solved analytically. It can however be solved numerically; we will return to this discussion when we later outline the necessary ingredients for the numerical solution of Eq. (62).
Keeping only the leading-twist terms in the expansion of the Wilson lines in Eq. (60) and performing the Gaussian averaging using Eq. (63), one can straightforwardly show that [8,29] φ qq,qq Similarly, and Substituting these leading twist distributions into Eq. (62), we find that the cross section simplifies greatly and can be expressed as The expression in Eq. (68) is the k ⊥ factorized expression for inclusive photon production in high energy QCD and is analogous to similar expressions derived previously for gluon production [40,7,22] and qq pair production [8] in this framework. To derive Eq. (69) we have taken advantage of the second line in the Eq. (46) to express T µ qq in terms of T µ q or T μ q . Alternately, it is useful to arrive at the results in Eqs. (68) and (69) starting directly from the perturbative computation of the amplitude. The leading twist diagrams are listed in Fig. 7. In fact, the leading twist amplitude in Lorenz gauge can immediately be read off from the original amplitudes (1)-(8) by expanding the Wilson lines to the first non-trivial order. In this case, it gives a non-zero contribution to the amplitude at the order O(ρ 1 p ρ 1 A ). The amplitudes (9)-(12) contain two insertions of the effective vertex (see Eq. (11)), and so they do not contribute at the order O(ρ 1 p ρ 1 A ). We can then write the leading twist amplitude as where the gluon fields are given by Fourier transforms of the expressions in Eqs. (3) and (7) and with the index β = 1, . . . , 8 corresponding to the respective diagram in Fig. 7 as Inserting the gluon fields from Eqs. (3) and (7), we can perform the light-cone integrations to obtain By taking the modulus squared and by performing the color averages as in Eqs. (57), (63) and (64), we can confirm that Eq. (68) is reproduced. The expression for the amplitude in Eq. (73), and in particular the individual contributions (72), can be compared to the amplitude for prompt photon hadroproduction first computed by Baranov et al. [9] and to the Z 0 hadroproduction from gluon fusion derived recently by Motyka et al. [10]. We have checked that our leading twist k ⊥ factorized result agrees exactly with Eq. (28) of Motyka et al.
One can also take the collinear limit, which as noted previously [41], can be obtained by taking k 1⊥ → 0 and k 2⊥ → 0 in the trace element, Θ(k 1⊥ , k 2⊥ )/k 2 1⊥ k 2 2⊥ , but not in the unintegrated functions. This means that the total external momenta vanishes, p ⊥ + q ⊥ + k γ⊥ = 0, thus guaranteeing momentum conservation. The limit is well defined thanks to the Ward identities ensuring that the amplitude vanishes linearly with k 1⊥ and k 2⊥ . The integration over k 1⊥ (k 2⊥ ) has to be taken then only over the unintegrated distribution functions, to get the inclusive photon cross section at fixed Q 2 and momentum fractions in the proton x p ∼ Q p S √ s e η kγ and in the nuclei where the |M gg→qqγ | 2 object is defined as An especially interesting point to note here is the sensitivity of the collinear result in Eq. (77) to the nuclear gluon distribution function f g,A . The CGC formalism used in this paper naturally extends beyond this result to the multi-parton distributions of the saturated nucleus as the momentum scales probed approach the semi-hard saturation scale Q A S .

Summary and Outlook
We have computed in this work the leading contribution to inclusive cross section for photon production at central rapidities in high energy proton-nucleus collisions. Our result for the triple differential cross section for a photon accompanied by a quark-antiquark pair is given in Eq. (61) and that for the inclusive photon cross section is given in Eq. (62). The result in this paper for the class III diagrams, in combination with the result for the class II diagrams in Fig. 2 obtained previously in [5], completes the NLO computation of inclusive photon production within the CGC framework.
The result in Eq. (62) has an identical structure to the expression for quark-antiquark pair production computed previously in [8]. Moreover as we show explicitly in Eq. (B.11) of Appendix B, the Low-Burnett-Kroll theorem tells us that, in the soft photon limit, the amplitude for photon production is precisely the amplitude for quark-antiquark pair production multiplied by a simple kinematic expression with the Lorentz structure of the photon.
The properties of the heavy quark pair production cross section derived in [8] have been explored extensively [28,29]. These results, when combined with a CGC+non-relativistic QCD (NRQCD) formalism [42], provide an excellent description of p+p J/Ψ production data [43] and p+A J/Ψ production data [44] at both RHIC and the LHC. The same formalism has also been employed to study J/Ψ production in p+A collisions within a color evaporation model of the hadronization of charm-anticharm pairs to J/Ψ mesons [45,46].
Computational techniques identical to those employed in these works can be used for quantitative estimates of inclusive photon production in p+A collisions at RHIC and LHC. The momentum dependent multi-parton correlations functions φ g,g A , φ qq,g A and φ qq,qq A appearing in Eq. (62) can be computed within the non-local Gaussian approximation we alluded to previously; this approximation reproduces a Langevin numerical implementation of the JIMWLK hierarchy [38] that delivers closed form results for such multi-parton correlators. The large N c limit of the multi-parton correlators provides an added simplification. Thus while the number of integrals to be performed are greater in the photon case relative to that of qq production, the computation of the former is feasible and will be reported on in future.
The numerical results will be very relevant for comparisons to anticipated results for direct photon measurements in p+A collisions at RHIC and the LHC. As noted, such measurements will be very sensitive to the nuclear gluon distribution function and will provide a quantitative estimate of power corrections to the same. Currently, the deuteron-gold measurements [30,31] at RHIC have covered the kinematic range 1 GeV k ⊥ 6 GeV with nearly real virtual photons and 5 GeV k ⊥ 16 GeV with real photons. It will be important to extend the latter measurements to lower k ⊥ to fully explore the gluon saturation regime. While the data is in agreement with pQCD, the remaining uncertainty allows for the possibility of thermal photons [47] on top of the pQCD results. Computations in our framework will help determine whether the latter are necessary.
Further, prompt photon measurements will provide an important test of the gluon saturation in general and, in particular, of the multi-parton correlators we have discussed here. They should corroborate the above mentioned studies of quarkonium production and may potentially be more robust since they are less sensitive to hadronization effects. We also note that our framework may provide insight into a long standing experimental puzzle [48] regarding the "anomalously" large photon production at soft momenta relative to predictions based on the Low-Burnett-Kroll theorem. We will investigate these issues in future work.
Another interesting feature of the results presented here is the comparative study of photon production in Lorentz gauge and light-cone gauge. While this study is a useful check of our results, it may also have further value in higher order computations that are feasible in this framework.
where β = 1, 2. To facilitate a comparison between different gauges we introduced extended tensors for Dirac indices as and for later convenience we also define Now we understand that we can express T µ g (k 1⊥ ) in the Lorenz gauge as Similarly, we proceed to evaluate the counterparts of regular diagrams (3)-(8) in Fig. 4. For example, let us look carefully at the evaluation of M µ 3 . The diagram (3) immediately leads to (A.10) in which only A µ p,LC contributes. In the above expression the k + 2 integration is trivial once we insert the effective vertex T (q + k γ , k 2 ), given in (11). The k − 2 integration picks up a singularity in S 0 (q + k γ − k 2 ). After the variable change change k 2 = P − k 1 we find We note that in the Lorenz gauge R µi 3 (k 1⊥ )(k 1i /P + ) would be replaced with −R µ− 3 . For diagrams (4)-(8) we can carry out similar procedures to define R µν β (k 1⊥ ) by generalizing γ − in R µ β (k 1⊥ ) to γ ν , and then M µ β is simply given by R µi β (k 1⊥ )(k 1i /P + ) instead of −R µ− β (k 1⊥ ). By analogy to Eq. (26) we define where T µ q (k 1⊥ ) and T μ q (k 1⊥ ) can be obtained as T µ q (k 1⊥ ) = R µ− q (k 1⊥ ) and T μ q (k 1⊥ ) = R µ− q (k 1⊥ ). Let us turn to contributions from diagrams (9)-(12) in Fig. 5. Here, we look specifically at the diagram (9) that yields, (A.14) The integrations over k + and k + 1 are trivial thanks to the nuclear effective vertex T . For the integrations over k − and k − 1 we pick up the singularities in the quark propagators S 0 (q +k γ −k) and in S 0 (q +k γ −k −k 1 ). We find .
(A. 23) In the first line we used (A.9) and in the second line we used the Dirac equation.
(A.27) This has exactly the same form as Eq. (A.23) but with an opposite sign. To end the proof we note that Eq. (A.27) is independent of k ⊥ and so we effectively have x ⊥ = y ⊥ in the fundamental Wilson lines in the amplitude (A.21). Then we can use the following identity, (A. 28) which will ensure that Eqs. (A.23) and (A.27) cancel out once these expressions are multiplied by the appropriate Wilson lines.

Appendix B. Properties of the photon production amplitude
In this Appendix we will prove that the amplitude satisfies the photon Ward identity and the famous Low-Burnett-Kroll soft photon theorem [13,14,15]. The Low-Burnett-Kroll theorem states that, in the infrared limit for the radiated photon momentum, i. e. k γ → 0, the amplitude should factorize into the product of the non-radiative amplitude, the photon polarization vector, and a vectorial structure that depends on momenta of the emitted charged particles.
We will rely on diagrammatic representations in the derivation in Sec. 2 though our conclusions will be valid of course for another derivation in Appendix A.
where the last piece defined by T g (k 1⊥ ) ≡ / C L (p + q, k 1⊥ ) (p + q) 2 , (B.8) is the expression obtained in the computation of the amplitude for quark-antiquark pair production [8].
For T µ qq , at the leading order, we only have to concern ourselves with diagrams (9) and (10) in Fig. 5. This is because the contributions from (11) and (12) do not contain any quark propagator in the form of Eq. (B.4), while (9) and (10) are divergent in the k γ → 0 limit, and this explains our identification of (9) and (10) as leading contributions and (11) and (12) as sub-leading ones. From this argument we can show that the T µ qq term with leading contributions only behaves as follows, * where the last piece defined by is the expression that also appears as a part of the amplitude for quark-antiquark pair production [8]. It should be noted that in the computation in Ref. [8] only the leading terms have been concerned, while the k γ → 0 limit of the amplitude in the present work can also pick up the sub-leading rest from diagrams (11) and (12). Combining both sets of leading contributions in the soft photon limit, we then get, * µ (k γ , λ)M µ (p, q, k γ ) → −q f e * µ (k γ , λ) where M(q, p) is the amplitude for qq pair production found in Ref. [8], and Eq. (B.11) verifies that the Low-Burnett-Kroll theorem in the k γ → 0 limit is satisfied.