Multiparticle production at mid-rapidity in the color-glass condensate

JHEP02(2019)024 Published for SISSA by Springer Received:September 26, 2018 Revised:December 6, 2018 Accepted:January 4, 2019 Published:February 5, 2019 Multiparticle production at mid-rapidity in the color-glass condensate Mauricio Martinez,aMatthew D. Sievertband Douglas E. Wertepnyc aDepartment of Physics, North Carolina State University, Raleigh, NC 27695, U.S.A. bTheoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, U.S.A. cDepartamento de F´ısica de Part´ıculas and IGFAE, Universidade de Santiago de Compostela, 15782 Santiago de Compostela, Galicia-Spain E-mail: [email protected],[email protected],[email protected] Abstract: In this paper, we compute a number of cross sections for the production of multiple particles at mid-rapidity in the semi-dilute / dense regime of the color-glass condensate (CGC) effective field theory. In particular, we present new results for the production of two quark-antiquark pairs (whether the same or different flavors) and for the production of one quark-antiquark pair and a gluon. We also demonstrate the existence of a simple mapping which transforms the cross section to produce a quark-antiquark pair into the corresponding cross section to produce a gluon, which we use to obtain various results and to cross-check them against the literature. We also discuss hadronization effects in the heavy flavor sector, writing explicit expressions for the production of various combinations of Dand ¯ Dmesons, J/ψ mesons, and light hadrons. The various multiparticle cross sections presented here contain a wealth of information and can be used to study heavy flavor production, charge-dependent correlations, and “collective” flow phenomena arising from initial-state dynamics. Keywords: Perturbative QCD, Quark-Gluon Plasma, Resummation ArXiv ePrint: 1808.04896 Open Access,c The Authors. Article funded by SCOAP3.https://doi.org/10.1007/JHEP02(2019)024 JHEP02(2019)024 Contents 1 Introduction 1 2 Production amplitudes for (anti)quark pairs and gluons 5 2.1 Quark / antiquark pair production amplitude 5 2.2 Gluon production amplitude 8 3 Cross sections for single pairs and gluons 10 3.1 Cross section for single-pair production 10 3.2 Cross section for single-gluon production 12 4 Double-inclusive cross sections for pairs and gluons 14 4.1 Cross section for double-pair production 14 4.1.1 Case 1: no fermion entanglement 15 4.1.2 Case 2: fermion entanglement 18 4.2 Cross section for quark + antiquark + gluon production 20 4.3 Double-gluon production 23 5 Hadronization in the heavy flavor sector 27 6 Conclusions 32 A Color averaging in the projectile and target 34 B Wilson line color algebra in momentum space 36 C Comparison with ref. [38]38 1 Introduction Correlations in the production of multiple soft or semi-hard particles in the mid-rapidity region of hadronic collisions are important probes of novel phenomena in quantum chromodynamics (QCD). Whether in proton-proton (pp), proton-nucleus (pA), or heavyion (AA) collisions, multiparticle production reflects the many-body correlations generated by QCD. In pp collisions, such correlations may be produced by quantum evolution through Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) evolution [1–3] or by a range of higher-order corrections to the hard part (see, e.g. [4,5]) or small-xevolution, including linear Balitsky-Fadin-Kuraev-Lipatov (BFKL) evolution [6,7], nonlinear Balitsky-Kovchegov (BK) evolution [8,9] and Jalilian-Marian-Iancu-McLerran-WeigertLeonidov-Kovner (JIMWLK) evolution [10–12]. The resulting correlations are sensitive probes of the perturbative hard vertex and of the strongly-ordered emission structure of the evolution equations. In pA collisions, these higher-order and evolution corrections are augmented by a new set of dynamical correlations arising from the enhancement of multiple scattering in the high charge densities of the heavy nucleus, characterized by the – 1 – JHEP02(2019)024 color-glass condensate (CGC) effective field theory (see e.g. [13] and references therein). The resulting correlations are sensitive probes of the multiple scattering dynamics, including the significant effects of Bose enhancement in the strong gluon fields [14,15]. Finally, in AA collisions (as well as potentially in high-multiplicity pp and pA collisions), all these initial state correlations are modified and complemented by the final-state dynamics of a strongly-coupled quark-gluon plasma (QGP) phase. A detailed characterization of multiparticle production in the strong color fields of the CGC is especially important in trying to differentiate intial-state effects from the finalstate dynamics of the QGP, where strongly-coupled interactions lead to substantial manybody correlations among soft and semi-hard particles. The canonical measures of this collective flow are the cumulants of azimuthal anisotropies vn{m}[16], with correlations among increasing numbers of particles reflected in higher values mof the cumulants. Other landmark properties believed to be possible in the QGP phase include the onset of novel transport mechanisms associated with the axial anomaly: the chiral magnetic effect, the chiral separation effect, the chiral vortical effect, and the chiral magnetic wave.1Signatures for all of these novel chiral dynamics are encoded in multiparticle correlations, often charge dependent, such as the same-sign and opposite-sign correlators γ112 and γ123 [18]. For all of these critical signatures of the quark-gluon plasma, it is essential to disentangle the “background” contributions coming from initial-state mechanisms to better quantify the properties of the QGP and improve the chances of discovering such novel anomalous dynamics. Accordingly, a substantial effort has been made in recent years to compute multiparticle production in the CGC framework. The purest realization of the CGC formalism is in the “dilute / dense” framework, in which density-enhanced effects of the “dilute projectile” are kept only to lowest order, while density-enhanced corrections in the “dense target” are resummed to all orders. Fewparticle production has been studied in the dilute / dense framework from the earliest days of the CGC formalism, starting with the inclusive single-gluon production cross section dσG[19–24] at mid-rapidity and followed shortly thereafter by the inclusive cross section dσq¯qof a single q¯qpair via gluon pair production [25–33]. Corrections to these production channels were also considered in the form of small-xevolution corrections [27,29]. However, a detailed computation of higher multiparticle production cross sections in the dilute / dense framework becomes increasingly difficult due to the proliferation of ways another soft particle could be radiated from a pre-existing one. A significant step toward overcoming this barrier was made through the development of the “semi-dilute / dense” framework [34]. This regime is designed to fill the gap between the dilute / dense regime, in which the projectile charge density is kept only to lowest order, and the dense / dense regime, where both projectile and target densities must be simultaneously resummed to all orders. The semi-dilute / dense framework is appropriate for “heavy-light ion collisions” intermediate to, say, pPb and PbPb collisions. For a collision between one light ion and one heavy ion, such as CuAu collisions, it is possible to construct a regime in which the large target density is resummed to all orders while corrections from 1For a review on chiral magnetic and vortical effects in high-energy nuclear collisions we refer to the reader to ref. [17]. – 2 – JHEP02(2019)024 the projectile density are calculated order by order in perturbation theory. In the semidilute / dense framework, higher-order corrections which are enhanced by the projectile density are more important than genuine quantum corrections. Formally, for a dense target nucleus with Anucleons and a semi-dilute projectile nucleus with anucleons, the semidilute / dense regime can be quantified by the hierarchy of scales α2 sA1/3∼ O(1) (1.1a) αsα2 sa1/31,(1.1b) or equivalently, in term of the saturation momenta Qs,p and Qs,t of the projectile and target respectively, Λ2 QCD Q2 s,p Q2 s,t.(1.2) With the help of the semi-dilute / dense framework, a number of significant steps have been taken in recent years toward the calculation of genuine multiparticle production in the CGC framework. The key simplification that makes this possible in the semi-dilute / dense framework is that the independent emission of new soft particles from the highdensity projectile becomes dominant over emission from the pre-existing system of soft particles. As such, the first observable computed in the semi-dilute / dense framework was the production cross section dσGG for two soft gluons [34,35]. A similar effort was made toward determining the production cross section dσqq for two quarks coming from separate q¯qpairs, with a partial calculation having been performed in ref. [36] emphasizing the new role played by Fermi-Dirac quantum statistics among the two pairs. This calculation later formed the basis of the partial calculation of the cross section dσqqG for two quark / antiquark pairs plus a gluon, with the intent of studying the CGC contribution to the same-sign correlators γ112 , γ123 [37]. It should be emphasized that in this important calculation [37], only correlations generated at the level of the wave functions were taken into account, without including the effects of multiple scattering that translate these wave functions into actual production cross sections. And very recently, a new attempt has been made to extend these calculations to the third order in the projectile charge density through the computation of the triple-gluon production cross section dσGGG [38]. Other notable developments in soft multiparticle production include the identification of Bose enhancement as a driving mechanism of the Ridge [14] in double-gluon production [14,15], the calculation of the soft double-photon cross section dσγγ [39], and the realization that soft double-pair production can be used to probe the gluon Wigner distribution with Weizs¨acker-Williams gauge structure [40]. A variant of the semi-dilute / dense power counting can be found in the form of the lowest-order “glasma graph” calculations, which have been used to calculate multiparticle correlations such as the triple-gluon cross section [41]. Other important developments in the calculation of multiparticle production in the CGC formalism have emphasized production in the forward regime, where the “hybrid factorization” framework makes it possible to rigorously relate the particle production cross sections to collinear parton distribution functions in the (semi-)dilute projectile, dressed with the effects of multiple scattering in the dense target [42–45]. In this approach, observables such as forward double valence-quark production cross sections dσqvqv[46], forward – 3 – JHEP02(2019)024 triple valence-quark production cross sections dσqvqvqv[47], and forward valence-quark + photon + gluon production [48] have been calculated. Similar studies of quadruple valencequark production cross sections dσqvqvqvqvhave also been considered in a “parton model” description [49,50] without the benefit of an underlying hybrid factorization. Other recent work on the subject also includes the demonstration [51] that two-gluon correlations can break the “accidental” back-to-back symmetry which occurs at lowest order and related phenomenology [52,53]. And finally, in a recent work [54], we have considered singleand double-pair production dσq¯qand dσ(q¯q)(q¯q)in coordinate space as a means of initializing spatial corrections of conserved charges in the quark-gluon plasma. In this paper, our primary goal is to systematically extend the calculation of multiparticle production at mid-rapidity in the semi-dilute / dense framework to higher orders. One of the key results we will derive here for the first time is the complete expression at exact Ncfor the double-pair production cross section dσ(q¯q)(q¯q)in momentum space, as written in eqs. (4.5), (4.11), and (4.16). This expression significantly generalizes the result obtained in ref. [36] by including contributions that were intentionally omitted there, by working with exact Nc, and by keeping the multiple scattering corrections to all orders. In a key conceptual development, we show in eq. (2.19) that it is possible to map the amplitude (and therefore, the cross section) for producing q¯qpairs into the corresponding quantities for producing gluons. Thus, we are able to directly map the double-pair cross section into the corresponding cross section dσ(q¯q)Gto produce a quark / antiquark pair and a gluon. This expression, as written in eq. (4.19), is also a new result. We also perform a number of validations of this gluonic mapping, verifying explicitly that it correctly reproduces the known results for singleand double-gluon production from the literature. While the preceding results all reflect the final-state production of multiple partons, they also open the door to a substantial program of computing hadronic-level observables derived from them. By convoluting the partonic-level results with the appropriate fragmentation functions or projection operators and long-distance matrix elements, we can translate these partonic-level cross sections to full hadronic cross sections. The resulting hadronic observables can be used to rigorously study the correlations among same-sign and opposite-sign charged hadrons, open and hidden heavy-flavor hadrons, heavy-flavor vs light hadrons, and more. The phenomenology based on these hadronic obserables will provide critical new insight into initial-state mechanisms for collective flow, quarkonium correlations, and charge-dependent correlations which form the background to anomalous chiral dynamics in the QGP. This paper is organized as follows. In section 2we construct the scattering amplitudes for the production of soft particles in momentum space, starting with the quark/antiquark production amplitude in section 2.1 and deriving the mapping to the gluon production amplitude in section 2.2. Then in section 3we compute the production cross section for a single q¯qpair in section 3.1 and map it in section 3.2 to the well-known gluon production cross section to validate the gluonic mapping. Then in section 4we proceed to calculate the new cross sections for the production of two sets of soft particles: double q¯qpair production in section 4.1, mixed q¯qG production in section 4.2, and double gluon production in section 4.3. The successful cross-check against the double-gluon production cross section – 4 – JHEP02(2019)024 k ¯ k=q−k q b u y x b σ σ′ Figure 1. The light-front wave functions to radiate a soft q¯qpair at mid-rapidity from a valence source, shown here as a quark. in section 4.3 represents another validation of the gluonic mapping derived in section 3.2. In section 5we utilize the techniques enumerated in ref. [55] to translate our partoniclevel cross sections into hadronic cross sections for the production of open and hidden heavy flavor as an illustration of how to straightforwardly apply the results derived here to hadronic observables. Finally, we conclude in section 6by reiterating the primary new results and exploring the many opportunities for phenomenological applications and further theoretical development which this work provides. In appendix Awe provide details of the Gaussian color averaging used for the (semi-)dilute projectile, in appendix Bwe formulate some useful algebraic properties of Wilson lines in momentum space, and in appendix Cwe point out the differences in the normalization of the double-gluon cross section against ref. [38]. Throughout this paper, we denote longitudinal momenta in light-front coordinates v±≡qg+− 2(v0±v3) and transverse vectors by v≡(v1 ⊥, v2 ⊥) with magnitudes vT≡ |v|. Different authors use different conventions for the light-front metric g+−; we will use g+−= 1, but it is also common to encounter g+−= 2. 2 Production amplitudes for (anti)quark pairs and gluons 2.1 Quark / antiquark pair production amplitude The amplitude to radiate a soft quark/antiquark pair at mid-rapidity has been derived many times in the literature [26,27]. In the notation of our previous work [54] as illustrated in figure 1, we denote the light-front wave functions [56,57] to radiate a soft q¯qpair as ψ1, ψ2, ψ3corresponding to the various time orderings of the scattering in the target fields. The term ψ1corresponds to scattering after the pair is created, ψ2to scattering after the gluon is emitted but before the pair is created, and ψ3to scattering before the pair is created. The three wave functions are not all independent, but satisfy ψ1+ψ2+ψ3= 0, and the explicit expressions are given by ψ1(q,κ) = −2gpα(1−α) κ2 T+m2+α(1−α)q2 T ×(δσ ,−σ02(1−α)+(1−2α)q·κ q2 T−iσ0q×κ q2 T−mσ0δσσ0"q1 ⊥ q2 T−iσ0q2 ⊥ q2 T#) (2.1a) – 5 – JHEP02(2019)024 ψ2(q,κ) = −2gpα(1−α) κ2 T+m2(δσ ,−σ0−(1−2α)q·κ q2 T +iσ0q×κ q2 T+mσ0δσσ0"q1 ⊥ q2 T−iσ0q2 ⊥ q2 T#) (2.1b) ψ3(q,κ) = −ψ1(q,κ)−ψ2(q,κ),(2.1c) where α≡k+ k++¯ k+is the fraction of the pair longitudinal momentum carried by the quark, qis the center-of-mass transverse momentum of the q¯qpair (i.e., the gluon), and κis the intrinsic transverse momentum of the quark splitting. In (2.1), we have omitted the explicit dependence of the wave functions on the momentum fraction αfor brevity. Note also that, in comparison to eqs. (21) of [54], we have removed a factor of the coupling gfrom the definition of the wave functions. This corresponds to absorbing this coupling constant into the scale µ2defined in (A.3) characterizing the sources of soft gluons. In terms of these wave functions, the single-pair amplitude summed over all time orderings is given in coordinate space by (see eqs. (30 - 31) of [54]) A(x, y, b)=(Vbta)hVxtaV† y−VbtaV† bψ1(u−b, x −y) +VutaV† u−VbtaV† bψ2(u−b, x −y)i,(2.2) where, as labeled in figure 1,x,y, and bare the final-state positions of the quark, antiquark, and valence quark, respectively, and u≡αx + (1 −α)yis the center-of-mass position of the q¯qpair (equal to the gluon position). The scattering of partons in the color fields of the target are described by Wilson lines in the fundamental or adjoint representations, Vx≡ Pexp ig Zdx+A−(x+,0−, x)(2.3a) Uab x≡Pexp ig Zdx+A− adj(x+,0−, x)ab ,(2.3b) where we work in the A+= 0 light cone gauge. With (2.2) written this way, the Wilson line Vbassociated with the valence quark will always cancel against a corresponding one in the complex-conjugate amplitude. It is convenient to translate the specific model of the projectile as a distribution of valence quarks into a generic continuous charge density. This can be accomplished by introducing the quantity ρa(b), which loosely corresponds to the wave function of a color source in the projectile at position bwhich radiates a soft gluon with color a. We can translate from the valence quark model of the projectile to the continuous color charge density by effectively replacing (Vbta)→ρa(b). (For another discussion of the translation between discrete and continuous charge distributions, see e.g. [51].) With this change of – 6 – JHEP02(2019)024 notation, we can Fourier transform the buildling block (2.2) into momentum space to obtain A(k, ¯ k) = Zd2x d2y d2b e−ik·xe−i¯ k·yρa(b) ×hVxtaV† y−VbtaV† bψ1(u−b, x −y) + VutaV† u−VbtaV† bψ2(u−b, x −y)i,(2.4) where the momenta of the final-state quark and antiquark are kand ¯ k, respectively. Note that the Fourier factor for the valence quark cancels because its position bis the same in the initial and final states under the eikonal approximation.2 Inserting the inverse transformation of the wave functions (2.1), Wilson lines, and source density ψi(u−b, x −y) = Zd2κ (2π)2 d2q (2π)2eiκ·(x−y)eiq·(u−b)ψi(q, κ) (2.5a) Vx=Zd2κ (2π)2eiκ·xV(κ) (2.5b) V† y=Zd2κ0 (2π)2e−iκ0·yV†(κ0) (2.5c) ρ(b) = Zd2` (2π)2ei`·bρ(`) (2.5d) gives A(k, ¯ k) = Zd2k0 (2π)2 d2q0 (2π)2ρa(q0)hV(k−k0)taV†(k−k0−q+q0)iψ1(q0,k0−αq0) +Zd2κ (2π)2 d2q0 (2π)2ρa(q0)hV(κ)taV†(κ−q+q0)iψ2(q0,k−αq) (2.6) −Zd2κ (2π)2 d2κ0 (2π)2ρa(q−κ+κ0)hV(κ)taV†(κ0)iψ1(q,k−αq)+ψ2(q,k−αq), with q=k+¯ kfor brevity. We can combine all three diagrams by redefining the dummy integration variables, obtaining the compact form A(k, ¯ k) = Zd2k0 (2π)2 d2¯ k0 (2π)2ρa(k0+¯ k0)hV(k−k0)taV†(¯ k0−¯ k)iΨ(k, ¯ k;k0,¯ k0),(2.7) where the differences among the three diagrams are all encoded in the combined wave function Ψ(k, ¯ k;k0,¯ k0)≡ψ1k0+¯ k0,(1 −α)k0−α¯ k0+ψ2k0+¯ k0,(1 −α)k−α¯ k −ψ1k+¯ k , (1 −α)k−α¯ k−ψ2k+¯ k , (1 −α)k−α¯ k.(2.8) 2At first glance, the amplitude (2.4) may appear to be problematic, because it contains an impact over impact parameters bof the source at the amplitude level, leading to two such impact parameter integrals in the cross section. This is true; however, when averaged over color states of the projectile as in (A.3), the correlator of two ρ’s possesses a delta function which sets these two positions equal. Thus the continuous charge distribution leads to one integral over d2bper source at the cross section level, as with the model of discrete valence quarks. – 7 – JHEP02(2019)024 q b x b Figure 2. The light-front wave function to radiate a soft gluon at mid-rapidity from a valence source, shown here as a quark. With this expression, it is easy to do the manipulations over all diagrams at once and particularly to study their color structure, since the Wilson lines enter in exactly the same form for all diagrams. As such, when we construct cross sections for the production of multiple pairs, we will only have to perform one calculation per diagrammatic topology, rather than having to repeat the calculation for many possible time orderings. These various topologies will correspond to different ways to contract the diagrams, including both the color matrix V taV†and the wave function Ψ, which is a matrix in the 2 ×2 spin space of the pair. 2.2 Gluon production amplitude In comparison with (2.2) for the production amplitude of a soft q¯qpair in coordinate space, the corresponding amplitude to emit a soft gluon is illustrated in figure 2and is given by Aa glue(x, b) = (Vbtb) (Ux)ab φ(x−b)−(taVb)φ(x−b),(2.9) where φis the light-front wave function φ(q)=2∗ λ·q q2 T δσvσ0 v(2.10a) φ(x−b) = i π ∗ λ·(x−b) (x−b)2 T δσvσ0 v(2.10b) to radiate a soft gluon from a valence quark projectile. Here σvand σ0 vare the spin states of the valence quark before and after gluon emission, and λis the spin of the emitted gluon. When we write the trace over the square of these wave functions, we mean the averaging over the quantum numbers of the initial state, together with a sum over the quantum numbers of the final state: trD[φ(q1)φ†(q2)] ≡1 2X λσvσ0 v φ(q1)φ∗(q2)=4q1·q2 q2 1Tq2 2T (2.11a) trD[φ(x)φ†(y)] ≡1 2X λσvσ0 v φ(x)φ∗(y) = 1 π2 x·y x2 Ty2 T .(2.11b) The first term of (2.9) corresponds to the shockwave passing through the gluon, the second term corresponds to the shockwave passing through the valence quark before the gluon is – 8 – JHEP02(2019)024 k2 ¯ k2 k′ 2 ¯ k′ 2 k′′ 2 ¯ k′′ 2 ρbρd k′ 2+¯ k′ 2k′′ 2+¯ k′′ 2 k1 ¯ k1 k′ 1 ¯ k′ 1 k′′ 1 ¯ k′′ 1 ρaρc k′ 1+¯ k′ 1k′′ 1+¯ k′′ 1 Figure 4. Double-pair production topologies without fermion entanglement, as calculated in eq. (4.6). dummy color indices a↔b. Squaring the full symmetrized amplitude (4.3) and converting to the cross section yields dσ(q¯q) (q¯q) d2k1dy1d2¯ k1d¯y1d2k2dy2d2¯ k2d¯y2 =1 2(2π)34 (4.5) ×DA(k1,¯ k1, k2,¯ k2) 2E−DA(k1,¯ k1, k2,¯ k2)A†(k1,¯ k2, k2,¯ k1)E+ (¯ k1↔¯ k2), where the notation +(¯ k1↔¯ k2) applies to both preceding terms. The task has now been reduced to calculating the two contributions in brackets: the case without fermion entanglement in the first term and the case with fermion entanglement in the second term. Both exercises are straightforward, and we calculate them in the following subsections. For maximum generality, we have considered here the case in which both produced pairs have the same flavor; if the flavor of the quark pairs is different, then all particles are distinguishable and only the first term of eq. (4.5) contributes. 4.1.1 Case 1: no fermion entanglement Squaring (4.4) for topologies with no fermion entanglement, as in figure 4, leads directly to DA(k1,¯ k1,k2,¯ k2) 2E=Z d 2{k0 1¯ k0 1k0 2¯ k0 2k00 1¯ k00 1k00 2¯ k00 2} ×Dρa(k0 1+¯ k0 1)ρb(k0 2+¯ k0 2)ρc∗(k00 1+¯ k00 1)ρd∗(k00 2+¯ k00 2)Eproj ×trDhΨ(k1,¯ k1;k0 1,¯ k0 1)Ψ†(k1,¯ k1;k00 1,¯ k00 1)itrDhΨ(k2,¯ k2;k0 2,¯ k0 2)Ψ†(k2,¯ k2;k00 2,¯ k00 2)i ×DtrchV(k1−k0 1)taV†(¯ k0 1−¯ k1)V(¯ k00 1−¯ k1)tcV†(k1−k00 1)i ×trchV(k2−k0 2)tbV†(¯ k0 2−¯ k2)V(¯ k00 2−¯ k2)tbV†(k2−k00 2)iEtgt.(4.6) – 15 – JHEP02(2019)024 There are now 3 possible “contractions” of the source colors, obtained in terms of (A.3): Dρa(k0 1+¯ k0 1)ρb(k0 2+¯ k0 2)ρc∗(k00 1+¯ k00 1)ρd∗(k00 2+¯ k00 2)Eproj = =δabδcd µ2k0 1+¯ k0 1+k0 2+¯ k0 2, k+ 1+¯ k+ 1+k+ 2+¯ k+ 2 ×µ2−k00 1−¯ k00 1−k00 2−¯ k00 2,−k+ 1−¯ k+ 1−k+ 2−¯ k+ 2 +δacδbd µ2k0 1+¯ k0 1−k00 1−¯ k00 1,0+µ2k0 2+¯ k0 2−k00 2−¯ k00 2,0+ +δadδbc µ2k0 1+¯ k0 1−k00 2−¯ k00 2, k+ 1+¯ k+ 1−k+ 2−¯ k+ 2 ×µ2k0 2+¯ k0 2−k00 1−¯ k00 1, k+ 2+¯ k+ 2−k+ 1−¯ k+ 1.(4.7) Note that only in the second term do the plus momenta combine to give 0+, even for this topology with no fermion entanglement. For a given set of color contractions, we will need to Fierz reduce the Wilson line traces of (4.6) twice. The algebra is straightforward, but it is convenient to define the following abbreviated notation: V1≡V(k1−k0 1)V5≡V(k2−k0 2) V† 2≡V†(¯ k0 1−¯ k1)V† 6≡V†(¯ k0 2−¯ k2) V3≡V(¯ k00 1−¯ k1)V7≡V(¯ k00 2−¯ k2) V† 4≡V†(k1−k00 1)V† 8≡V†(k2−k00 2).(4.8) With this shorthand, the Wilson line tensor entering (4.6) is Ωabcd 1= trchV1taV† 2V3tcV† 4itr hV5tbV† 6V7tdV† 8i,(4.9) and we can straightforwardly compute the various contraction of (4.7): δabδcd Ωabcd 1=N2 c 4Dˆ D4(1674) ˆ D4(5238)E−1 4D8(16785234) −1 4D8(12385674) + 1 4Dˆ D4(1234) ˆ D4(5678)E,(4.10a) δacδbd Ωabcd 1=N4 c 4Dˆ D2(32) ˆ D2(14) ˆ D2(76) ˆ D2(58)E−N2 c 4Dˆ D2(32) ˆ D2(14) ˆ D4(5678)E −N2 c 4Dˆ D4(1234) ˆ D2(76) ˆ D2(58)E+1 4Dˆ D4(1234) ˆ D4(5678)E,(4.10b) δadδbc Ωabcd 1=N2 c 4Dˆ D4(1854) ˆ D4(7236)E−1 4D8(34185672) −1 4D8(12367854) + 1 4Dˆ D4(1234) ˆ D4(5678)E,(4.10c) with the numbers in parentheses denoting the arguments of the corresponding Wilson lines in (4.8). These various traces are ilustrated in figure 5. Fantastically, we only have to do one calculation per topology (i.e., source contraction) because all of the different time orderings enter on the same footing. Combining all these terms back into (4.6) yields the – 16 – JHEP02(2019)024 h ˆ D2 ˆ D2 ˆ D2 ˆ D2i h ˆ D2 ˆ D2 ˆ D4i h ˆ D2 ˆ D6i h ˆ D4 ˆ D4iD8 Figure 5. Illustration of the Wilson line traces hˆ D2ˆ D2ˆ D2ˆ D2i,hˆ D2ˆ D2ˆ D4i,hˆ D2ˆ D6i,hˆ D4ˆ D4i, and D8contributing to double-pair production in eqs. (4.11) and (4.16). complete result for topologies with no fermion entanglement: DA(k1,¯ k1,k2,¯ k2) 2E=N4 c 4Z d 2{k0 1¯ k0 1k0 2¯ k0 2k00 1¯ k00 1k00 2¯ k00 2} ×trDhΨ(k1,¯ k1;k0 1,¯ k0 1)Ψ†(k1,¯ k1;k00 1,¯ k00 1)itrDhΨ(k2,¯ k2;k0 2,¯ k0 2)Ψ†(k2,¯ k2;k00 2,¯ k00 2)i ×µ2(k0 1+¯ k0 1+k0 2+¯ k0 2, k+ 1+¯ k+ 1+k+ 2+¯ k+ 2)µ2(−k00 1−¯ k00 1−k00 2−¯ k00 2,−k+ 1−¯ k+ 1−k+ 2−¯ k+ 2) ×1 N2 cDˆ D4(k1−k0 1,¯ k0 2−¯ k2,¯ k00 2−¯ k2,k1−k00 1)ˆ D4(k2−k0 2,¯ k0 1−¯ k1,¯ k00 1−¯ k1,k2−k00 2)E −1 N4 cD8(k1−k0 1,¯ k0 2−¯ k2,¯ k00 2−¯ k2,k2−k00 2,k2−k0 2,¯ k0 1−¯ k1,¯ k00 1−¯ k1,k1−k00 1) −1 N4 cD8(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1,k2−k00 2,k2−k0 2,¯ k0 2−¯ k2,¯ k00 2−¯ k2,k1−k00 1) +1 N4 cDˆ D4(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1,k1−k00 1)ˆ D4(k2−k0 2,¯ k0 2−¯ k2,¯ k00 2−¯ k2,k2−k00 2)E +µ2(k0 1+¯ k0 1−k00 1−¯ k00 1,0+)µ2(k0 2+¯ k0 2−k00 2−¯ k00 2,0+) ×Dˆ D2(¯ k00 1−¯ k1,¯ k0 1−¯ k1)ˆ D2(k1−k0 1,k1−k00 1)ˆ D2(¯ k00 2−¯ k2,¯ k0 2−¯ k2)ˆ D2(k2−k0 2,k2−k00 2)E −1 N2 cDˆ D2(¯ k00 1−¯ k1,¯ k0 1−¯ k1)ˆ D2(k1−k0 1,k1−k00 1)ˆ D4(k2−k0 2,¯ k0 2−¯ k2,¯ k00 2−¯ k2,k2−k00 2)E −1 N2 cDˆ D4(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1,k1−k00 1)ˆ D2(¯ k00 2−¯ k2,¯ k0 2−¯ k2)ˆ D2(k2−k0 2,k2−k00 2)E +1 N4 cDˆ D4(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1,k1−k00 1)ˆ D4(k2−k0 2,¯ k0 2−¯ k2,¯ k00 2−¯ k2,k2−k00 2)E +µ2(k0 1+¯ k0 1−k00 2−¯ k00 2, k+ 1+¯ k+ 1−k+ 2−¯ k+ 2)µ2(k0 2+¯ k0 2−k00 1−¯ k00 1, k+ 2+¯ k+ 2−k+ 1−¯ k+ 1) ×1 N2 cDˆ D4(k1−k0 1,k2−k00 2,k2−k0 2,k1−k00 1)ˆ D4(¯ k00 2−¯ k2,¯ k0 1−¯ k1,¯ k00 1−¯ k1,¯ k0 2−¯ k2)E −1 N4 cD8(¯ k00 1−¯ k1,k1−k00 1,k1−k0 1,k2−k00 2,k2−k0 2,¯ k0 2−¯ k2,¯ k00 2−¯ k2,¯ k0 1−¯ k1) −1 N4 cD8(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1,¯ k0 2−¯ k2,¯ k00 2−¯ k2,k2−k00 2,k2−k0 2,k1−k00 1) +1 N4 cDˆ D4(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1,k1−k00 1)ˆ D4(k2−k0 2,¯ k0 2−¯ k2,¯ k00 2−¯ k2,k2−k00 2)E (4.11) – 17 – JHEP02(2019)024 k1 k2 ¯ k2 ¯ k′′ 1 k′′ 2 ρc k′′ 2+¯ k′′ 1 ¯ k1 k′ 2 ¯ k′ 2 k′ 2+¯ k′ 2 ρb k′ 1 ¯ k′ 1 k′ 1+¯ k′ 1 ρa k′′ 1 ¯ k′′ 2 k′′ 1+¯ k′′ 2 ρd Figure 6. Double-pair production “Pac Man” topologies with fermion entanglement, as calculated in eq. (4.12). 4.1.2 Case 2: fermion entanglement The interference term in (4.5) contains the topologies with the fermions being entangled, such that the pairs in the amplitude “swap ownership” of one of the fermions in going to the complex-conjugate amplitude. This contribution arises from the Fermi-Dirac statistics of identical particles and therefore does not contribute if the pairs have different flavors. Taking the interference of (4.4) as shown in the “Pac Man” type diagram of figure 6we have directly DA(k1,¯ k1, k2,¯ k2)A†(k1,¯ k2, k2,¯ k1)E=Z d 2{k0 1¯ k0 1k0 2¯ k0 2k00 1¯ k00 1k00 2¯ k00 2} ×Dρa(k0 1+¯ k0 1)ρb(k0 2+¯ k0 2)ρc∗(k00 2+¯ k00 1)ρd∗(k00 1+¯ k00 2)Eproj ×trDhΨ(k1,¯ k1;k0 1,¯ k0 1) Ψ†(k2,¯ k1;k00 2,¯ k00 1) Ψ(k2,¯ k2;k0 2,¯ k0 2)Ψ†(k1,¯ k2;k00 1,¯ k00 2)i ×DtrchV(k1−k0 1)taV†(¯ k0 1−¯ k1)V(¯ k00 1−¯ k1)tcV†(k2−k00 2) ×V(k2−k0 2)tbV†(¯ k0 2−¯ k2)V(¯ k00 2−¯ k2)tdV†(k1−k00 1)iEtgt.(4.12) In the same way as (4.7), we form the 3 contractions of the projectile sources: Dρa(k0 1+¯ k0 1)ρb(k0 2+¯ k0 2)ρc∗(k00 2+¯ k00 1)ρd∗(k00 1+¯ k00 2)Eproj = (4.13) =δabδcd µ2(k0 1+¯ k0 1+k0 2+¯ k0 2, k+ 1+¯ k+ 1+k+ 2+¯ k+ 2) ×µ2(−k00 2−¯ k00 1−k00 1−¯ k00 2,−k+ 2−¯ k+ 1−k+ 1−¯ k+ 2) +δacδbdµ2(k0 1+¯ k0 1−k00 2−¯ k00 1, k+ 1−k+ 2)µ2(k0 2+¯ k0 2−k00 1−¯ k00 2, k+ 2−k+ 1) +δadδbcµ2(k0 1+¯ k0 1−k00 1−¯ k00 2,¯ k+ 1−¯ k+ 2)µ2(k0 2+¯ k0 2−k00 2−¯ k00 1,¯ k+ 2−¯ k+ 1). – 18 – JHEP02(2019)024 Using the same shorthand notation as (4.8), the Wilson line tensor entering (4.12) is Ωabcd 2= trc[V1taV† 2V3tcV† 8V5tbV† 6V7tdV† 4],(4.14) and we can compute the various color contractions in the same way: δabδcd Ωabcd 2=Nc 4D8(16785234) −Nc 4Dˆ D4(3852) ˆ D4(1674)E −Nc 4Dˆ D4(5678) ˆ D4(1234)E+1 4NcD8(12385674).(4.15a) δacδbd Ωabcd 2=N3 c 4Dˆ D2(32) ˆ D2(76) ˆ D4(1854)E−Nc 4Dˆ D2(32) ˆ D6(185674)E −Nc 4Dˆ D2(76) ˆ D6(123854)E+1 4NcD8(12385674).(4.15b) δadδbc Ωabcd 2=N3 c 4Dˆ D2(58) ˆ D2(14) ˆ D4(7236)E−Nc 4Dˆ D2(14) ˆ D6(385672)E −Nc 4Dˆ D2(58) ˆ D6(123674)E+1 4NcD8(12385674).(4.15c) Combining these back into (4.12) yields the complete result for topologies with fermion entanglement: DA(k1,¯ k1,k2,¯ k2)A†(k1,¯ k2,k2,¯ k1)E=N3 c 4Z d 2{k0 1¯ k0 1k0 2¯ k0 2k00 1¯ k00 1k00 2¯ k00 2} ×trDhΨ(k1,¯ k1;k0 1,¯ k0 1)Ψ†(k2,¯ k1;k00 2,¯ k00 1)Ψ(k2,¯ k2;k0 2,¯ k0 2)Ψ†(k1,¯ k2;k00 1,¯ k00 2)i ×µ2(k0 1+¯ k0 1+k0 2+¯ k0 2, k+ 1+¯ k+ 1+k+ 2+¯ k+ 2)µ2(−k00 2−¯ k00 1−k00 1−¯ k00 2,−k+ 2−¯ k+ 1−k+ 1−¯ k+ 2) ×h1 N2 cD8(k1−k0 1,¯ k0 2−¯ k2,¯ k00 2−¯ k2,k2−k00 2,k2−k0 2,¯ k0 1−¯ k1,¯ k00 1−¯ k1,k1−k00 1) −1 N2 cDˆ D4(¯ k00 1−¯ k1,k2−k00 2,k2−k0 2,¯ k0 1−¯ k1)ˆ D4(k1−k0 1,¯ k0 2−¯ k2,¯ k00 2−¯ k2,k1−k00 1)E −1 N2 cDˆ D4(k2−k0 2,¯ k0 2−¯ k2,¯ k00 2−¯ k2,k2−k00 2)ˆ D4(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1,k1−k00 1)E +1 N4 cD8(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1,k2−k00 2,k2−k0 2,¯ k0 2−¯ k2,¯ k00 2−¯ k2,k1−k00 1)i +µ2(k0 1+¯ k0 1−k00 2−¯ k00 1, k+ 1−k+ 2)µ2(k0 2+¯ k0 2−k00 1−¯ k00 2, k+ 2−k+ 1) ×hDˆ D2(¯ k00 1−¯ k1,¯ k0 1−¯ k1)ˆ D2(¯ k00 2−¯ k2,¯ k0 2−¯ k2)ˆ D4(k1−k0 1,k2−k00 2,k2−k0 2,k1−k00 1)E −1 N2 cDˆ D2(¯ k00 1−¯ k1,¯ k0 1−¯ k1)ˆ D6(k1−k0 1,k2−k00 2,k2−k0 2,¯ k0 2−¯ k2,¯ k00 2−¯ k2,k1−k00 1)E −1 N2 cDˆ D2(¯ k00 2−¯ k2,¯ k0 2−¯ k2)ˆ D6(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1,k2−k00 2,k2−k0 2,k1−k00 1)E +1 N4 cD8(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1,k2−k00 2,k2−k0 2,¯ k0 2−¯ k2,¯ k00 2−¯ k2,k1−k00 1)i +µ2(k0 1+¯ k0 1−k00 1−¯ k00 2,¯ k+ 1−¯ k+ 2)µ2(k0 2+¯ k0 2−k00 2−¯ k00 1,¯ k+ 2−¯ k+ 1) ×hDˆ D2(k2−k0 2,k2−k00 2)ˆ D2(k1−k0 1,k1−k00 1)ˆ D4(¯ k00 2−¯ k2,¯ k0 1−¯ k1,¯ k00 1−¯ k1,¯ k0 2−¯ k2)E −1 N2 cDˆ D2(k1−k0 1,k1−k00 1)ˆ D6(¯ k00 1−¯ k1,k2−k00 2,k2−k0 2,¯ k0 2−¯ k2,¯ k00 2−¯ k2,¯ k0 1−¯ k1)E −1 N2 cDˆ D2(k2−k0 2,k2−k00 2)ˆ D6(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1,¯ k0 2−¯ k2,¯ k00 2−¯ k2,k1−k00 1)E +1 N4 cD8(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1,k2−k00 2,k2−k0 2,¯ k0 2−¯ k2,¯ k00 2−¯ k2,k1−k00 1)i(4.16) – 19 – JHEP02(2019)024 The expression for the cross section (4.5), together with the two classes of topologies (4.11) and (4.16) constitute the first complete and exact solution to the 4-particle inclusive (q¯q) (q¯q) cross section at this order. These expressions are one of the primary results of this paper. 4.2 Cross section for quark + antiquark + gluon production In the same way as in section 3.2, we can now take the “gluonic limit” of the q¯qpair k2,¯ k2→1 2q2in the double-pair expression (4.11). Note that, once we replace one of the pairs with a gluon, there are no longer any possible fermion entanglement topologies, since all three final state particles q¯qG are now distinguishable. Thus we need only take the limit of the topologies in (4.11), immediately obtaining: dσ(q¯q)G d2k1dy1d2¯ k1d¯y1d2q2dy2=1 2(2π)33N4 c 4Z d 2{k0 1¯ k0 1k00 1¯ k00 1q0 2q00 2δk0 2δk00 2} ×trDhΨ(k1,¯ k1;k0 1,¯ k0 1)Ψ†(k1,¯ k1;k00 1,¯ k00 1)itrDhΦ(q2;q0 2)Φ†(q2;q00 2)i ×µ2(k0 1+¯ k0 1+q0 2, k+ 1+¯ k+ 1+q+ 2)µ2(−k00 1−¯ k00 1−q00 2,−k+ 1−¯ k+ 1−q+ 2) ×1 N2 cDˆ D4(k1−k0 1,1 2q0 2−1 2q2−δk0 2,1 2q00 2−1 2q2−δk00 2,k1−k00 1) ׈ D4(1 2q2−1 2q0 2−δk0 2,¯ k0 1−¯ k1,¯ k00 1−¯ k1,1 2q2−1 2q00 2−δk00 2)E −1 N4 cD8(k1−k0 1,1 2q0 2−1 2q2−δk0 2,1 2q00 2−1 2q2−δk00 2,1 2q2−1 2q00 2−δk00 2, 1 2q2−1 2q0 2−δk0 2,¯ k0 1−¯ k1,¯ k00 1−¯ k1,k1−k00 1) −1 N4 cD8(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1,1 2q2−1 2q00 2−δk00 2, 1 2q2−1 2q0 2−δk0 2,1 2q0 2−1 2q2−δk0 2,1 2q00 2−1 2q2−δk00 2,k1−k00 1) +1 N4 cDˆ D4(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1,k1−k00 1) ׈ D4(1 2q2−1 2q0 2−δk0 2,1 2q0 2−1 2q2−δk0 2,1 2q00 2−1 2q2−δk00 2,1 2q2−1 2q00 2−δk00 2)E +µ2(k0 1+¯ k0 1−k00 1−¯ k00 1,0+)µ2(q0 2−q00 2,0+) ×Dˆ D2(¯ k00 1−¯ k1,¯ k0 1−¯ k1)ˆ D2(k1−k0 1,k1−k00 1)ˆ D2(1 2q00 2−1 2q2−δk00 2,1 2q0 2−1 2q−δk0 2) ׈ D2(1 2q2−1 2q0 2−δk0 2,1 2q2−1 2q00 2−δk00 2)E −1 N2 cDˆ D2(¯ k00 1−¯ k1,¯ k0 1−¯ k1)ˆ D2(k1−k0 1,k1−k00 1) ׈ D4(1 2q2−1 2q0 2−δk0 2,1 2q0 2−1 2q2−δk0 2,1 2q00 2−1 2q2−δk00 2,1 2q2−1 2q00 2−δk00 2)E −1 N2 cDˆ D4(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1,k1−k00 1)ˆ D2(1 2q00 2−1 2q2−δk00 2,1 2q0 2−1 2q2−δk0 2) ׈ D2(1 2q2−1 2q0 2−δk0 2,1 2q2−1 2q00 2−δk00 2)E +1 N4 cDˆ D4(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1,k1−k00 1) ׈ D4(1 2q2−1 2q0 2−δk0 2,1 2q0 2−1 2q2−δk0 2,1 2q00 2−1 2q2−δk00 2,1 2q2−1 2q00 2−δk00 2)E – 20 – JHEP02(2019)024 +µ2(k0 1+¯ k0 1−q00 2, k+ 1+¯ k+ 1−q+ 2)µ2(q0 2−k00 1−¯ k00 1, q+ 2−k+ 1−¯ k+ 1) ×1 N2 cDˆ D4(k1−k0 1,1 2q2−1 2q00 2−δk00 2,1 2q2−1 2q0 2−δk0 2,k1−k00 1) ׈ D4(1 2q00 2−1 2q2−δk00 2,¯ k0 1−¯ k1,¯ k00 1−¯ k1,1 2q0 2−1 2q2−δk0 2)E −1 N4 cD8(¯ k00 1−¯ k1,k1−k00 1,k1−k0 1,1 2q2−1 2q00 2−δk00 2, 1 2q2−1 2q0 2−δk0 2,1 2q0 2−1 2q2−δk0 2,1 2q00 2−1 2q2−δk00 2,¯ k0 1−¯ k1) −1 N4 cD8(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1,1 2q0 2−1 2q2−δk0 2, 1 2q00 2−1 2q2−δk00 2,1 2q2−1 2q00 2−δk00 2,1 2q2−1 2q0 2−δk0 2,k1−k00 1) +1 N4 cDˆ D4(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1,k1−k00 1) ׈ D4(1 2q2−1 2q0 2−δk0 2,1 2q0 2−1 2q2−δk0 2,1 2q00 2−1 2q2−δk00 2,1 2q2−1 2q00 2−δk00 2)E, (4.17) where as before we have changed variables into q0 2≡k0 2+¯ k0 2and δk0 2≡1 2(k0 2−¯ k0 2), and similarly for q00 2and δk00 2. After integration over δk0 2, δk00 2, most of these terms vanish. These integrals act only on the interaction terms, and as we saw in eqs. (3.8) and (3.9), any time the same momentum δk(0,00) 2appears in adjacent arguments of the same trace, those Wilson lines will cancel. In canceling, they lead to delta functions of the momentum difference following (B.5), which is always either δ2(q0 2−q2) or δ2(q00 2−q2), and these terms drop out due to the vanishing of the wave function as in (3.9). The only interaction terms which do not vanish correspond to lines 1, 5, 7, and 9 out of the 12 terms in braces. As in (3.10), these terms instead simplify the Wilson line traces by causing some of the coordinate-space arguments to be repeated, with the composite object being Fourier transformed to momentum space. One of the two surviving operators was already calculated in (3.10): the square of the dipole amplitude. The other nonvanishing operator is a partial reduction of the double quadrupole (see figure 7): Z d 2{δk0 2δk00 2}Dˆ D4(p1,1 2q0 2−1 2q2−δk0 2,1 2q00 2−1 2q2−δk00 2, p2) ׈ D4(1 2q2−1 2q0 2−δk0 2, p3, p4,1 2q2−1 2q00 2−δk00 2)E =Zd2{x1y1x2y2z1w1z2w2}Z d 2{δk0 2δk00 2} ×e−ip1·x1ei1 2q0 2−1 2q2−δk0 2·y1e−i1 2q00 2−1 2q2−δk00 2·x2eip2·y2 ×e−i1 2q2−1 2q0 2−δk0 2·z1eip3·w1e−ip4·z2ei1 2q2−1 2q00 2−δk00 2·w2 ×Dˆ D4(x1, y1, x2, y2)ˆ D4(z1, w1, z2, w2)E =Zd2{x1y1x2y2w1z2}e−ip1·x1ei(q0 2−q2)·y1e−i(q00 2−q2)·x2eip2·y2eip3·w1e−ip4·z2 ×Dˆ D4(x1, y1, x2, y2)ˆ D4(y1, w1, z2, x2)E ≡D4,4(p1, q0 2−q2, q00 2−q2, p2;p3p4).(4.18) – 21 – JHEP02(2019)024 h ˆ D2 ˆ D2| ˆ D2|2iD4,4 h ˆ D4| ˆ D2|2i Figure 7. Illustration of the Wilson line traces hˆ D2ˆ D2|ˆ D2|2i,hˆ D4|ˆ D2|2iand D4,4contributing to q¯qG production in eq. (4.19). The result is the complete expression for the q¯qG cross section in momentum space, with the corresponding operators illustrated in figure 7: dσ(q¯q)G d2k1dy1d2¯ k1d¯y1d2q2dy2=1 2(2π)33N4 c 4Z d 2{k0 1¯ k0 1k00 1¯ k00 1q0 2q00 2} ×trDhΨ(k1,¯ k1;k0 1,¯ k0 1) Ψ†(k1,¯ k1;k00 1,¯ k00 1)itrDhΦ(q2;q0 2) Φ†(q2;q00 2)i ×µ2(k0 1+¯ k0 1+q0 2, k+ 1+¯ k+ 1+q+ 2)µ2(−k00 1−¯ k00 1−q00 2,−k+ 1−¯ k+ 1−q+ 2) ×h1 N2 cD4,4(k1−k0 1, q0 2−q2, q00 2−q2, k1−k00 1;¯ k0 1−¯ k1¯ k00 1−¯ k1)i +µ2(k0 1+¯ k0 1−k00 1−¯ k00 1,0+)µ2(q0 2−q00 2,0+) ×Dˆ D2(¯ k00 1−¯ k1,¯ k0 1−¯ k1)ˆ D2(k1−k0 1, k1−k00 1) ˆ D2 2(q00 2−q2, q0 2−q2)E −1 N2 cDˆ D4(k1−k0 1,¯ k0 1−¯ k1,¯ k00 1−¯ k1, k1−k00 1) ˆ D2 2(q00 2−q2, q0 2−q2)E +µ2(k0 1+¯ k0 1−q00 2, k+ 1+¯ k+ 1−q+ 2)µ2(q0 2−k00 1−¯ k00 1, q+ 2−k+ 1−¯ k+ 1) ×1 N2 cD4,4(k1−k0 1, q2−q00 2, q2−q0 2, k1−k00 1;¯ k0 1−¯ k1¯ k00 1−¯ k1).(4.19) To our knowledge, eq. (4.19) represents the first complete calculation of q¯qG production in the CGC framework, and this compact form in momentum space comprises an exact solution at this order, at finite Nc. We emphasize, however, that this cross section applies only for the production of a quark and antiquark of the same flavor. This new expression is the second primary result of this paper. – 22 – JHEP02(2019)024 4.3 Double-gluon production As a final cross-check of the preceding calculations, let us use the mapping (2.19) on the remaining q¯qpair in (4.19) to obtain the double-gluon production cross section, which we can compare with explicit results in the literature. As before, we take the limit k1,¯ k1→1 2q1 in (4.19), adjust the count of 1 2(2π)3in the prefactor, and change integration variables to q0 1≡k0 1+¯ k0 1and δk0 1≡1 2(k0 1−¯ k0 1), and similarly for q00 1and δk00 1. This gives dσGG d2q1dy1d2q2dy2=1 2(2π)32N4 c 4Z d 2{q0 1q00 1q0 2q00 2δk0 1δk00 1}(4.20) ×trDhΦ(q1;q0 1)Φ†(q1;q00 1)itrDhΦ(q2;q0 2)Φ†(q2;q00 2)i ×µ2(q0 1+q0 2, q+ 1+q+ 2)µ2(−q00 1−q00 2,−q+ 1−q+ 2) ×h1 N2 cD4,4(1 2q1−1 2q0 1−δk0 1,q0 2−q2,q00 2−q2,1 2q1−1 2q00 1−δk00 1;1 2q0 1−1 2q1−δk0 11 2q00 1−1 2q1−δk00 1)i +µ2(q0 1−q00 1,0+)µ2(q0 2−q00 2,0+) ×Dˆ D2(1 2q00 1−1 2q1−δk00 1,1 2q0 1−1 2q1−δk0 1)ˆ D2(1 2q1−1 2q0 1−δk0 1,1 2q1−1 2q00 1−δk00 1) × ˆ D2 2(q00 2−q2,q0 2−q2)E −1 N2 cDˆ D4(1 2q1−1 2q0 1−δk0 1,1 2q0 1−1 2q1−δk0 1,1 2q00 1−1 2q1−δk00 1,1 2q1−1 2q00 1−δk00 1) × ˆ D2 2(q00 2−q2,q0 2−q2)E +µ2(q0 1−q00 2, q+ 1−q+ 2)µ2(q0 2−q00 1, q+ 2−q+ 1) ×h1 N2 cD4,4(1 2q1−1 2q0 1−δk0 1,q2−q00 2,q2−q0 2,1 2q1−1 2q00 1−δk00 1;1 2q0 1−1 2q1−δk0 11 2q00 1−1 2q−δk00 1)i. Of the four interaction terms remaining in the braces, the third one (quadrupole trace) vanishes after integration over δk0 1δk00 1due to repeated adjacent arguments; the second one (double dipole) is the same as (3.10); and the first and last ones are further reductions of the double-quadrupole (4.18): Z d 2{δk0 1δk00 1}D4,4(1 2q1−1 2q0 1−δk0 1,p1,p2,1 2q1−1 2q00 1−δk00 1;1 2q0 1−1 2q1−δk0 11 2q00 1−1 2q1−δk00 1) =Zd2{x1y1x2y2w1z2}Z d 2{δk0 1δk00 1}e−i1 2q1−1 2q0 1−δk0 1·x1eip1·y1e−ip2·x2ei1 2q1−1 2q00 1−δk00 1·y2 ×ei1 2q0 1−1 2q1−δk0 1·w1e−i1 2q00 1−1 2q1−δk00 1·z2Dˆ D4(x1,y1,x2,y2)ˆ D4(y1,w1,z2,x2)E =Zd2{x1y1x2y2w1z2}e−i(q1−q0 1)·x1eip1·y1e−ip2·x2ei(q1−q00 1)·y2 ˆ D4(x1,y1,x2,y2) 2 ≡|D4|2(q1−q0 1,p1,p2,q1−q00 1),(4.21) – 23 – JHEP02(2019)024 which is just the square of the quadrupole. Thus the double-gluon cross section takes the especially compact form in momentum space dσGG d2q1dy1d2q2dy2=1 2(2π)32N4 c 4Z d 2{q0 1q00 1q0 2q00 2}(4.22) ×trDhΦ(q1;q0 1)Φ†(q1;q00 1)itrDhΦ(q2;q0 2)Φ†(q2;q00 2)i ×µ2(q0 1+q0 2, q+ 1+q+ 2)µ2(−q00 1−q00 2,−q+ 1−q+ 2)h1 N2 c|D4|2(q1−q0 1,q0 2−q2,q00 2−q2,q1−q00 1)i +µ2(q0 1−q00 1,0+)µ2(q0 2−q00 2,0+)D ˆ D2 2(q00 1−q1,q0 1−q1) ˆ D2 2(q00 2−q2,q0 2−q2)E +µ2(q0 1−q00 2, q+ 1−q+ 2)µ2(q0 2−q00 1, q+ 2−q+ 1)h1 N2 c|D4|2(q1−q0 1,q2−q00 2,q2−q0 2,q1−q00 1)i). At this point, we can directly compare the cross section (4.22) against the known expressions in the literature; ref. [38] gives the GG cross section in momentum space, and ref. [34] gives the cross section in coordinate space. Here we will perform the transformation to coordinate space to demonstrate exact agreement with ref. [34], and in appendix Cwe present the additional comparison with ref. [38] in momentum space. Most of the work in performing the cross-check against ref. [34] comes from unfolding the compact expression (4.22) back into coordinate space. Inserting the Fourier transforms of the various quantities gives dσGG d2q1dy1d2q2dy2 =1 2(2π)32N4 c 4Zd2{x1y1x2y2x0 1y0 1x0 2y0 2b1b2}(4.23) ×trDhΦ(x1−b1;x0 1−b1)Φ†(y1−b1;y0 1−b1)itrDhΦ(x2−b2;x0 2−b2)Φ†(y2−b2;y0 2−b2)i ×D ˆ D2 2(x0 1,y0 1) ˆ D2 2(x0 2,y0 2)E ×e−iq1·(x0 1−y0 1+x1−y1)e−iq2·(x0 2−y0 2+x2−y2)µ2(b1,0+)µ2(b2,0+) +trDhΦ(x1−b1;x0 1−b1)Φ†(y2−b2;y0 2−b2)itrDhΦ(x2−b2;x0 2−b2)Φ†(y1−b1;y0 1−b1)i ×1 N2 c|D4|2(x0 1,y0 1,x0 2,y0 2) ×e−iq1·(x0 1−y0 2+x1−y2−b1+b2)e−iq2·(x0 2−y0 1+x2−y1−b2+b1)µ2(b1, q+ 1−q+ 2)µ2(b2, q+ 2−q+ 1) +e−iq1·(x0 1−y0 2+x1−y2−b1+b2)eiq2·(x0 2−y0 1+x2−y1+b1−b2)µ2(b1, q+ 1+q+ 2)µ2(b2,−q+ 1−q+ 2). In arriving at (4.23), we have used symmetry properties of the squared dipole and the gluon wave functions (2.11). Inserting (3.13) in the four wave functions generates 16 terms, each containing 4 delta functions. As before, after integrating over those delta functions, we drop any remaining primes on the integration variables; this leads to all of the various wave functions in a given set of terms having the same arguments, with the differences residing in the Wilson line interactions. – 24 – JHEP02(2019)024 the intrinsic transverse momentum of the pojectile is characterized by the projectile saturation scale Qs,p and is computed order by order in perturbation theory; and the intrinsic transverse momentum of hadronization is simply of order ΛQCD and is not enhanced by any high density scales. Thus for the semi-dilute / dense regime, we have the hierarchy of scales (1.2), for which we can neglect the intrinsic transverse momentum characterized by the fragmentation functions. This hierarchy of scales is in the same spirit as the framework of hybrid factorization, in which the production of two distinguishable heavy quarkonia has recently been calculated [40]. One potential drawback to the asymmetric treatment of the projectile, target, and fragmentation sectors is that the fragmentation functions employed above do not possess an unambiguous factorization scale µF. This feature also applies to the description of fragmentation employed in ref. [55]. While the particular fragmentation functions given in ref. [61] do not explicitly refer to a factorization scale µF, fragmentation functions in general — and the available fragmentation functions for light hadrons in particular — carry such a dependence on an arbitrary scale. This arbitrary scale dependence in the fragmentation sector is compensated by the scale dependence of the parton distribution functions of the projectile and target, such that the observable cross section is overall invariant under the renormalization-group evolution of its various nonperturbative pieces. In our case, the projectile, target, and fragmentation sectors have all been treated very differently, such that the scale dependence coming from the projectile and target distributions is ambiguous.4In a more complete treatment which puts the nonperturbative projectile, target, and fragmentation sectors on comparable footing, the cancellation of this factorization scale dependence would become explicit. This is what is seen explicitly in the standard hybrid factorization framework of the dilute / dense regime [64,65], and we expect that the same would be true for the semi-dilute / dense regime in a treatment such as that of ref. [40], which formulates hybrid factorization in terms of double parton distributions of the projectile. Finally, let us note that while the ICEM is a particularly simple and convenient model for describing the hadronization of a quarkonium state from a q¯qpair, it is by no means unique. A variety of other descriptions of this hadronization process exist, in particular the effective field theory of Non-Relativistic Quantum Chromodynamics (NRQCD) [66]. While different hadronization formalisms have their own advantages and disadvantages, NRQCD has the particular advantage of being a self-consistent effective field theory of QCD. Employing it requires a more detailed treatment than just a simple convolution of the partonic c¯ccross section as done above, using a series of projection operators to select out the quantum numbers of the c¯cstate appropriate for a given hadronization channel. Of these projection operators, the color projections onto singlet and octet c¯cstates will require a more detailed implementation because they will modify the Wilson lines which enter the multipole traces. These various projections are in principle straightforward, but they are beyond the scope of this paper; we leave the incorporation of an NRQCD-based approach 4This is not to be confused with the scale dependence on a rapidity regulator; the RG evolution in rapidity is contained within the rapidity dependence of the Wilson line traces. This dependence is generally characterized by the JIMWLK functional evolution equation. – 31 – JHEP02(2019)024 to quarkonium production for future work. It is interesting to note, however, that inclusive J/ψ production at small xin the CGC framework was studied in ref. [55] using both NRQCD and the ICEM for hadronization, concluding that the ICEM is a good approximation to the NRQCD approach due to the dominance of the 3S[8] 1production channel. 6 Conclusions In this paper, we have computed a number of cross sections for the production of multiple particles at mid-rapidity in the semi-dilute / dense regime of the color-glass condensate effective field theory. At the partonic level, we have computed for the first time the production cross sections for two quark/antiquark pairs (q¯q) (q¯q) (eqs. (4.5), (4.11), and (4.16)) and for one quark/antiquark pair plus a gluon (q¯q)G(eq. (4.19)). The double-pair expression significantly extends previous work [36] in that it is fully differential in all four particles, includes all time orderings and all orders of multiple rescattering in the target fields, and is evaluated with exact Nc. These new partonic cross sections are one of the primary new results of this paper. Additionally, we proved a simple mapping (2.19) between the production amplitude for a q¯qpair and the production amplitude for a gluon, which we used to obtain the (q¯q)Gcross section from the (q¯q) (q¯q) cross section, and which we validated by crosschecking the single-gluon (3.11) and double-gluon (4.22) production cross sections against the literature [13,34,38]. The mapping (2.19) and its application to derive whole classes of cross sections from the multi-pair cross section is the second primary result of this paper. Finally, in section 5we discussed how to translate the partonic cross sections computed here into hadronic ones in the heavy flavor sector through the use of collinear fragmentation functions for open heavy flavor and the Improved Color Evaporation Model [63] for heavy quarkonia. This procedure translates each of the partonic cross sections into a range of hadronic observables, allowing us to write down expressions for the production of: (D¯ D) (D¯ D) — eq. (5.2); (D¯ D) (J/ψ) — eq. (5.6); (J/ψ) (J/ψ) — eq. (5.5); (D¯ D)h— eq. (5.3); and (J/ψ)h— eq. (5.7), where his any light hadron. These expressions open the door for a wide range of phenomenology to study the production and correlations of many particles in the heavy flavor sector, and they are the third primary result of this paper. The ability to perform a small number of ab initio calculations in the CGC formalism at the partonic level to simultaneously predict the production cross sections and correlations of a wide variety of hadronic observables has the potential to broadly test the initial-state mechanisms as an explanation for the observed correlations in heavyand heavy-light ion collisions. Genuine multiparticle production computed within the CGC framework makes it possible to self-consistently compute higher cumulants such as v2{4}and charge-dependent correlations like γ112 from purely initial-state mechanisms. Similarly, the coordinate-space program begun in ref. [54] aspires to take partonic correlations such as these as inputs to the initial conditions of subsequent hydrodynamic evolution, including contributions from conserved charges due to quark production. As we continue to extend this program of genuine multiparticle production, beyond simple approximations like the so-called “dilute / dilute glasma graphs” or the large-Ncapproximation, we anticipate that it will open up – 32 – JHEP02(2019)024 broad opportunities to test the effects of initial-state correlations, with and without the impact of strongly-coupled final-state interactions. One potential barrier to such a comprehensive program of multi-hadron phenomenology in the CGC approach is that multiparticle cross sections, like the ones calculated here, invoke higher and higher n-point correlators of Wilson lines, up to the octupole D8for double-pair production. These operators become increasingly difficult to evaluate, even in simple models like the MV model, for which analytic expressions are only available for the 4-point functions [67]. However, a compelling argument summarized in ref. [38] and attributed to ref. [47] suggests that, up to corrections suppressed by the large area of the target, n-point correlators can in general be factorized into products of 2-point correlators (dipoles), which are well-constrained in theory and phenomenology. If this argument holds, then the increasingly complex Wilson line structure is no obstacle to the pursuit of multihadron phenomenology. Aside from the applications already discussed above, there are a number of other direct extensions of this method we can pursue in future work. One is to repeat the double-pair calculation of section 4.1 in coordinate space to study the spatial correlations among quarks and antiquarks; as discussed in ref. [54], these double-pair correlations are the dominant effect for same-sign charged particles and for opposite-sign charged particles at distances larger than the inverse quark mass 1/m. Another is to extend the calculations performed here for double-pair production to triple-pair production: (q¯q) (q¯q) (q¯q). The number of permutations will increase substantially in going to triple-pair production, but the fundamental mechanics of the calculation will not change, and the compact expression (2.7) in momentum space makes such an extension tractable. Moreover, the gluonic mapping (2.19) will make it possible to immediately translate the triple-pair cross section into a whole family of related cross sections:(q¯q) (q¯q) (q¯q); (q¯q) (q¯q)G; (q¯q)G G; and G G G. Last, we note that the expressions for hadronization in the heavy flavor sector we explore in section 5can be significantly improved and extended. The heavy flavor sector is convenient as a justification for the assumption that a given hadron (like a Dmeson) in the final state is dominated by fragmentation from a given parton (like a cquark). In principle, a sum over all partonic channels with appropriate fragmentation functions will relax this assumption and make it possible to study fragmentation into several identified hadrons. As we extend the program to compute multiparticle production in the CGC framework, we will include more and more of these partonic channels, allowing a complete calculation of correlations in inclusive hadron production. In particular, the CGC contribution to the charge-dependent correlations γ112 and γ123 which form the background to the signal for the chiral magnetic effect is of special importance. While some exploratory work on this subject was done in ref. [37], it includes neither scattering in the target fields nor hadronization, both of which are likely to strongly modify the charge-dependent correlations. However, with the calculation of a range of partonic channels and appropriate charge-dependent fragmentation functions, a robust computation of the hadronic charge correlations becomes possible. For all of these reasons, we believe that the theoretical advancements presented in this paper represent a significant step toward implementing a comprehensive program of phenomenology to study multi-hadron correlations from the initial state. – 33 – JHEP02(2019)024 Acknowledgments The authors wish to thank N. Armesto, D. Pitonyak, J. Noronha-Hostler, V. Skokov, P. Tribedy, and R. Venugopalan for useful discussions. This work is supported in part by the U.S. Department of Energy grant DE-FG02-03ER41260 and the BEST (Beam Energy Scan Theory) DOE Topical Collaboration (MM), DOE Contract No. DE-AC5206NA25396 and the DOE Early Career Program (MS), the European Research Council 39 grant HotLHC ERC-2011-StG-279579, Ministerio de Ciencia e Innovaci´on of Spain under project FPA2014-58293-C2-1-P and Unidad de Excelencia Mar´ıa de Maeztu under project MDM-2016-0692, Xunta de Galicia (Conseller´ıa de Educaci´on) and FEDER (DW). A Color averaging in the projectile and target In calculating cross sections, we will need to compute squares and interferences of the elementary building block (2.7) and average them over various fluctuating quantities. Aside from the averaging over the quantum numbers of the produced particles, which is performed in the usual way, the event averaging covers three distinct types of fluctuations. These are fluctuations in the color fields of the projectile, fluctuations in the color fields of the target, and fluctuations over the global collision geometry. These three types of averaging factorize from one another, such that we can average the sources ρin the color fields of the projectile, which we denote h···iproj, separately from the averaging of the Wilson lines over the color fields in the target, which we denote h···itgt or simply h···iwhen there is no ambiguity. The collision geometry is characterized by an impact parameter Bbetween the centers of the projectile and target, which can be either held fixed at the cross section level or integrated out at the end of the calculation. Similarly, there may be other parameters describing the overall collision geometry, such as the angular orientation of a non-spherical nucleus like uranium; these global parameters will also be integrated out at the cross section level. In this paper, we make no particular assumption about the nature of the averaging in the color fields of the target; our final expressions for the cross sections will involve various traces of Wilson lines (2.3) corresponding to color dipole, quadrupole, sextupole, and octupole operators. We denote those corresponding operators by ˆ D2(x, y)≡1 Nc tr hVxV† yi(A.1a) ˆ D4(x, y, z, w)≡1 Nc tr hVxV† yVzV† wi(A.1b) ˆ D6(x, y, z, w, u, v)≡1 Nc tr hVxV† yVzV† wVuV† vi(A.1c) ˆ D8(x, y, z, w, u, v, r, s)≡1 Nc tr hVxV† yVzV† wVuV† vVrV† si.(A.1d) For quantities that have already been averaged we drop the hat over the operator, writing e.g. D2(x, y) = hˆ D2(x, y)i. When computing similar traces with adjoint Wilson lines, we will denote the operator with the superscript “adj”, writing e.g. ˆ Dadj 2(x, y)≡1 N2 c−1Uab xU† yba . – 34 – JHEP02(2019)024 And we will maintain the same notation whether invoking Wilson lines in coordinate space or momentum space, writing e.g. ˆ D2(p, q) = 1 Nctr[V(p)V†(q)] in momentum space. It is also important to note that the (averaged) Wilson line traces from eqs. (A.1) implicitly depend on a rapidity scale Y. This rapidity scale regulates the light-cone divergences associated with higher-order corrections to these operators, and there are a range of different schemes available to regulate these divergences. Physically, we can think of this scale as being set by the total rapidity interval Y∝ln s ⊥2of the collision, and the quantum evolution with the running of this scale is given by the JIMWLK evolution equations [10–12] (or the large-Ncanalogue, the BK equation [8,9]). This evolution is triggered when the rapidity interval is parametrically large, Y∼1 αs. On the other hand, we restrict ourselves here to the case when the produced particles are close enough in rapidity that we do not need to consider quantum evolution in the rapidity interval between the particles. Formally, this means ∆yij <1 αsfor the rapidity interval ∆yij between any two particles i, j tagged at mid-rapidity. For the average h···iproj over projectile color fields, we’ll use a Gaussian averaging procedure inspired by the McLerran-Venugopalan (MV) model [68]. Gaussian averaging corresponds to limiting the interaction of a source particle to two gluons, which is the lowest order in perturbation theory that preserves color neutrality. For Gaussian color averaging, the expectation value of products of several ρ’s factorizes into a sum over all possible pairwise “contractions,” such that it is only necessary to specify the two-point function hρ ρito fully specify the result of the averaging procedure: hρa(x)ρb(y)ρc∗(z)ρd∗(w)iproj =Dρa(x)ρb(y)Eproj Dρc∗(z)ρd∗(w)Eproj (A.2) +hρa(x)ρc∗(z)iproj Dρb(y)ρd∗(w)Eproj +Dρa(x)ρd∗(w)Eproj Dρb(y)ρc∗(z)Eproj . The original MV model [68] has been generalized in a number of different ways in the literature; for our purposes, it is useful to enumerate three distinct physical assumptions about the two-point function: Dρa(x)ρb∗(y)Eproj =       δab δ2(x−y)δ(x−−y−)µ2(x, x−) Locality δab δ(x−−y−)µ2(|x−y|2 T, x−) 2D Transl.Inv. δab δ2(x−y)δ(x−−y−)µ2(x−) Both (A.3a) Dρa(q1)ρb∗(q2)Eproj =       δab µ2(q1−q2, q+ 1−q+ 2) Locality δab (2π)2δ2(q1−q2)µ2(q2 1T, q+ 1−q+ 2) 2D Transl.Inv. δab (2π)2δ2(q1−q2)µ2(q+ 1−q+ 2) Both. (A.3b) For averages with both of the color fields in the amplitude hρ ρi, the momentum of the second source is reversed: q2→ −q2, and for averages with both of the color fields in the complex conjugate amplitude hρ∗ρ∗i, the momentum of the first source is reversed: q1→ −q1. The first case in eqs. (A.3), “Locality,” is the one we will employ throughout this paper. This assumption describes color fluctuations characterized by Gaussian random noise: they – 35 – JHEP02(2019)024 are only correlated locally at the same point, with different points in space having totally uncorrelated random fluctuations. The second case, “2D Translational Invariance,” does not require locality, but does assume that the average distribution of source charges is uniform in the transverse plane; this assumption is employed in references such as [14] and [36]. The simplest case uses “Both” locality and 2D translational invariance; this assumption is the one employed by the original MV model [68] and generalized somewhat in ref. [69]. An interesting argument about the general structure of the two-point function requiring only very weak assumptions about color neutrality was recently given in ref. [46]; for additional discussion about the properties of the two-point function see e.g. refs. [69] and [38]. Whatever physical assumptions are made about the two-point function in the Gaussian averaging, the color charge density fluctuations are characterized by the scale µ2, which is related to the saturation scale of the projectile Q2 s,proj ∝µ2[13,70]. As defined in (A.3), the scale µ2effectively contains a factor of the coupling g2associated with radiating a soft gluon from the color sources; some references prefer to write this factor explicitly, but we will use the conventions of (A.3) in which that coupling is contained within µ2. Note also that the dimensions of µ2as written in (A.3) are different among the different physical assumptions. It should also be emphasized that the averaging performed here for the projectile, denoted h···iproj, refers only to averaging over color configurations at a fixed collision geometry. The scale µ2in (A.3) is written with a dependence on the transverse position x, which implicitly keeps fixed the global impact parameter Bbetween the projectile and target. In translating back from the continuous color charge densities used in (A.3) to the discrete valence quark distributions used for example in ref. [54], the corresponding dictionary is µ2(b1)···µ2(bn)→Zd2Bg2 2Nc Tproj(b1−B)···g2 2Nc Tproj(bn−B),(A.4) where the factor of g2accounts for the coupling constant in the emission of the soft gluon from the valence quark source and the 1 2Ncarises from averaging over the color states of the valence quark: 1 Nctrc[tatb] = 1 2Ncδab. B Wilson line color algebra in momentum space The cancellation of Wilson lines which is trivial in coordinate space takes on a more subtle form in momentum space: 1 = VxV† x=Zd2p (2π)2 d2q (2π)2ei(p−q)·xV(p)V†(q) (B.1) Clearly we can’t just cancel momentum-space Wilson lines of equal argument on the righthand side: V(p)V†(p)6= 1. In fact, there are two actions taking place on the momentumspace Wilson lines resulting in the cancellation: a Fourier transform over the relative momentum (p−q) and integral over the center-of-mass momentum p+q 2. The only way – 36 – JHEP02(2019)024 for these two integrals to lead to unity on the left-hand is if one of these actions results in a delta function, which the other one picks up. But at this level it is not clear which operation should generate the delta function. On the other hand, we can engineer a cancellation of Wilson lines in momentum space by explicitly generating Wilson lines in coordinate space with the same argument: V(p)V†(p+q) = Zd2x d2y e−ip·xei(p+q)·yVxV† y.(B.2) Clearly if we integrate over the shared momentum p, we generate a delta function that sets x=yand cancels the Wilson lines. The remaining integral over d2ythen generates a new delta function: Zd2p (2π)2V(p)V†(p+q) = Zd2y eiq·yVyV† y= (2π)2δ2(q).(B.3) This condition is necessary, and as it turns out, it is also sufficient to guarantee (B.1). Using (B.3) in (B.1) we obtain VxV† x=Zd2p (2π)2 d2q (2π)2e−iq·xV(p)V†(p+q) =Zd2q (2π)2e−iq·x(2π)2δ2(q) = 1,(B.4) so we can conclude that the two conditions are in fact equivalent: Zd2p (2π)2V(p)V†(p+q) = (2π)2δ2(q)↔hVxV† x= 1i.(B.5) Another important feature is the conversion between Wilson line traces in the fundamental and adjoint representations. In coordinate space, the squares of fundamental dipoles and quadrupoles can be directly converted into adjoint dipoles and quadrupoles, plus a constant term which is Ncsuppressed: ˆ Dadj 2(x,y)≡1 N2 c−1Uab xU† yba =Nc 2CF ˆ D2(x,y) 2−1 N2 c−1(B.6a) ˆ Dadj 4(x,y,z,w)≡1 N2 c−1Uab xU† ybc Ucd zU† wda =Nc 2CF ˆ D4(x,y,z,w) 2−1 N2 c−1.(B.6b) In going to momentum space, that additive constant becomes instead proportional to delta functions: ˆ Dadj 2(p, q) = Nc 2CF ˆ D2 2(p, q)−1 N2 c−1(2π)4δ2(p)δ2(q) (B.7a) ˆ Dadj 4(p, q, p0, q0) = Nc 2CF ˆ D4 2(p, q, p0, q0)−1 N2 c−1(2π)8δ2(p)δ2(q)δ2(p0)δ2(q0).(B.7b) – 37 – JHEP02(2019)024 C Comparison with ref. [38] Given the successful cross-check in section 4.3 of our double-gluon cross section (4.22) against ref. [34] in coordinate space, it is also desirable to compare against the expressions given in eqs. (11 - 13) of [38] in momentum space. There, they keep maximally general µ2=µ2(q1, q2) for the projectile, the Lipatov vertex in their eq. (14) corresponds to our wave functions (2.18) and (2.11) tr[Φ(q;q0) Φ†(q;q00)] = 8Li ⊥(q, q −q0)Li ⊥(q, q −q00),(C.1) and they use the expressions for the adjoint dipole and quadrupole traces from (B.7). When used in (4.22), the delta function terms of (B.7) always vanish, because they lead to a vanishing wave function Φ(q;q) = 0. Using this in (4.22) gives dσ d2q1dy1d2q2dy2 =16 [2(2π)3]2Z d 2{q0 1q00 1q0 2q00 2} ×Li ⊥(q1, q1−q0 1)Li ⊥(q1, q1−q00 1)Lj ⊥(q2, q2−q0 2)Lj ⊥(q2, q2−q00 2) ×((N2 c−1) µ2(q0 1+q0 2)µ2(−q00 1−q00 2)Dadj 4(q1−q0 1, q0 2−q2, q00 2−q2, q1−q00 1) + (N2 c−1)2µ2(q0 1−q00 1)µ2(q0 2−q00 2)Dˆ Dadj 2(q00 1−q1, q0 1−q1)ˆ Dadj 2(q00 2−q2, q0 2−q2)E + (N2 c−1) µ2(q0 1−q00 2)µ2(q0 2−q00 1)Dadj 4(q1−q0 1, q2−q00 2, q2−q0 2, q1−q00 1)),(C.2) which can be compared with their eqs. (11-13). The three lines of (C.2) are referred to as A, B, and C, and are given in eqs. (11), (12), and (13) of [38], respectively. A direct comparison is obtained by changing variables to correspond to the arguments of the Wilson line traces, and the final result is dσ d2q1dy1d2q2dy2 = (4παs)216 [2(2π)3]2Z d 2{`1`2`3`4} ×((N2 c−1)Dadj 4(`1,`2,`3,`4)µ2(q1+q2−`1+`4) g2 µ2(−q1−q2−`3+`4) g2 ×Li ⊥(q1,`1)Li ⊥(q1,`2)Lj ⊥(q2,−`3)Lj ⊥(q2,−`4) +(N2 c−1)2Dˆ Dadj 2(`1,`2)ˆ Dadj 2(`3,`4)Eµ2(`2−`1) g2 µ2(`4−`3) g2 ×Li ⊥(q1,`1)Li ⊥(q1,`2)Lj ⊥(q2,`3)Lj ⊥(q2,`4) +(N2 c−1)Dadj 4(`1,`2,`3,`4)µ2(q1−q2−`1+`4) g2 µ2(q2−q1+`2−`3) g2 ×Li ⊥(q1,`1)Li ⊥(q1,`2)Lj ⊥(q2,`3)Lj ⊥(q2,`4)),(C.3) where we have used various symmetry properties of the momentum-space adjoint traces. Eq. (C.3) agrees perfectly with eqs. (11-13) of ref. [38], except for the prefactor. The – 38 – JHEP02(2019)024 conversion factor (4παs g2)2is just a trivial difference of convention: we take our µ2to include the gluon-emission coupling g2, while they insert the coupling explicitly. We do, however, differ by the numerically significant prefactor 16 [2(2π)3]2≈6.5×10−5which arises from the invariant phase space of the two-gluon final state. While the precise value of the prefactor is largely unimportant for the physical content of the calculation, we note that including this prefactor is necessary to obtain absolute agreement with ref. [34]. Numerically, inserting this prefactor suppresses the cross section by more than 4 orders of magnitude, which will surely be important for the phenomenology of multiparticle production. Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited. References [1] Y.L. Dokshitzer, Calculation of the Structure Functions for Deep Inelastic Scattering and e+e−Annihilation by Perturbation Theory in Quantum Chromodynamics.,Sov. Phys. JETP 46 (1977) 641 [INSPIRE]. [2] V.N. Gribov and L.N. Lipatov, Deep inelastic e p scattering in perturbation theory,Sov. J. Nucl. Phys. 15 (1972) 438 [INSPIRE]. [3] G. Altarelli and G. Parisi, Asymptotic Freedom in Parton Language,Nucl. Phys. B 126 (1977) 298 [INSPIRE]. [4] R.K. Ellis, H. Georgi, M. Machacek, H.D. Politzer and G.G. Ross, Perturbation Theory and the Parton Model in QCD,Nucl. Phys. B 152 (1979) 285 [INSPIRE]. [5] G. Curci, W. Furmanski and R. Petronzio, Evolution of Parton Densities Beyond Leading Order: The Nonsinglet Case,Nucl. Phys. B 175 (1980) 27 [INSPIRE]. [6] E.A. Kuraev, L.N. Lipatov and V.S. Fadin, The Pomeranchuk Singularity in Nonabelian Gauge Theories,Sov. Phys. JETP 45 (1977) 199 [INSPIRE]. [7] I.I. Balitsky and L.N. Lipatov, The Pomeranchuk Singularity in Quantum Chromodynamics, Sov. J. Nucl. Phys. 28 (1978) 822 [INSPIRE]. [8] Y.V. Kovchegov, Small-x F2structure function of a nucleus including multiple Pomeron exchanges,Phys. Rev. D 60 (1999) 034008 [hep-ph/9901281] [INSPIRE]. [9] I. Balitsky, Operator expansion for high-energy scattering,Nucl. Phys. B 463 (1996) 99 [hep-ph/9509348] [INSPIRE]. [10] J. Jalilian-Marian, A. Kovner, L.D. McLerran and H. Weigert, The intrinsic glue distribution at very small x,Phys. Rev. D 55 (1997) 5414 [hep-ph/9606337] [INSPIRE]. [11] J. Jalilian-Marian, A. Kovner, A. Leonidov and H. Weigert, The Wilson renormalization group for low x physics: Towards the high density regime,Phys. Rev. D 59 (1998) 014014 [hep-ph/9706377] [INSPIRE]. [12] E. Iancu, A. Leonidov and L.D. McLerran, The renormalization group equation for the color glass condensate,Phys. Lett. B 510 (2001) 133 [hep-ph/0102009] [INSPIRE]. [13] Y.V. Kovchegov and E. Levin, Quantum Chromodynamics at High Energy, Cambridge University Press, (2012). – 39 – JHEP02(2019)024 [14] T. Altinoluk, N. Armesto, G. Beuf, A. Kovner and M. Lublinsky, Bose enhancement and the ridge,Phys. Lett. B 751 (2015) 448 [arXiv:1503.07126] [INSPIRE]. [15] T. Altinoluk, N. Armesto and D.E. Wertepny, Correlations and the ridge in the Color Glass Condensate beyond the glasma graph approximation,JHEP 05 (2018) 207 [arXiv:1804.02910] [INSPIRE]. [16] M. Luzum and H. Petersen, Initial State Fluctuations and Final State Correlations in Relativistic Heavy-Ion Collisions,J. Phys. G 41 (2014) 063102 [arXiv:1312.5503] [INSPIRE]. [17] D.E. Kharzeev, J. Liao, S.A. Voloshin and G. Wang, Chiral magnetic and vortical effects in high-energy nuclear collisions — A status report,Prog. Part. Nucl. Phys. 88 (2016) 1 [arXiv:1511.04050] [INSPIRE]. [18] CMS collaboration, Constraints on the chiral magnetic effect using charge-dependent azimuthal correlations in pPb and PbPb collisions at the CERN Large Hadron Collider, Phys. Rev. C 97 (2018) 044912 [arXiv:1708.01602] [INSPIRE]. [19] Y.V. Kovchegov and A.H. Mueller, Gluon production in current-nucleus and nucleon-nucleus collisions in a quasiclassical approximation,Nucl. Phys. B 529 (1998) 451 [hep-ph/9802440] [INSPIRE]. [20] J.P. Blaizot, F. Gelis and R. Venugopalan, High-energy pA collisions in the color glass condensate approach. 1. Gluon production and the Cronin effect,Nucl. Phys. A 743 (2004) 13 [hep-ph/0402256] [INSPIRE]. [21] Y.V. Kovchegov and K. Tuchin, Inclusive gluon production in DIS at high parton density, Phys. Rev. D 65 (2002) 074026 [hep-ph/0111362] [INSPIRE]. [22] A. Dumitru and L.D. McLerran, How protons shatter colored glass,Nucl. Phys. A 700 (2002) 492 [hep-ph/0105268] [INSPIRE]. [23] B.Z. Kopeliovich, A.V. Tarasov and A. Schafer, Bremsstrahlung of a quark propagating through a nucleus,Phys. Rev. C 59 (1999) 1609 [hep-ph/9808378] [INSPIRE]. [24] J. Jalilian-Marian and Y.V. Kovchegov, Saturation physics and deuteron-Gold collisions at RHIC,Prog. Part. Nucl. Phys. 56 (2006) 104 [hep-ph/0505052] [INSPIRE]. [25] E.M. Levin, M.G. Ryskin, Yu. M. Shabelski and A.G. Shuvaev, Heavy quark production in semihard nucleon interactions,Sov. J. Nucl. Phys. 53 (1991) 657 [INSPIRE]. [26] J.P. Blaizot, F. Gelis and R. Venugopalan, High-energy pA collisions in the color glass condensate approach. 2. Quark production,Nucl. Phys. A 743 (2004) 57 [hep-ph/0402257] [INSPIRE]. [27] Y.V. Kovchegov and K. Tuchin, Production of q¯qpairs in proton-nucleus collisions at high energies,Phys. Rev. D 74 (2006) 054014 [hep-ph/0603055] [INSPIRE]. [28] F. Gelis and R. Venugopalan, Large mass q¯qproduction from the color glass condensate, Phys. Rev. D 69 (2004) 014019 [hep-ph/0310090] [INSPIRE]. [29] H. Fujii, F. Gelis and R. Venugopalan, Quark pair production in high energy pA collisions: General features,Nucl. Phys. A 780 (2006) 146 [hep-ph/0603099] [INSPIRE]. [30] F. Gelis, K. Kajantie and T. Lappi, Chemical thermalization in relativistic heavy ion collisions,Phys. Rev. Lett. 96 (2006) 032304 [hep-ph/0508229] [INSPIRE]. – 40 –