Multiparticle Production at Mid-Rapidity in the Color-Glass Condensate

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 $D$ and $\bar D$ mesons, $J/\psi$ 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.


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][2][3] or by a range of higher-order corrections to the hard part (see, e.g. [4,5]) or small-x evolution, including linear Balitsky-Fadin-Kuraev-Lipatov (BFKL) evolution [6,7], nonlinear Balitsky-Kovchegov (BK) evolution [8,9] and Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) evolution [10][11][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 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 v n {m} [16], with correlations among increasing numbers of particles reflected in higher values m of 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. 1 Signatures 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][20][21][22][23][24] at mid-rapidity and followed shortly thereafter by the inclusive cross section dσ qq of a single qq pair via gluon pair production [25][26][27][28][29][30][31][32][33]. Corrections to these production channels were also considered in the form of small-x evolution 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 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 A nucleons and a semi-dilute projectile nucleus with a nucleons, the semidilute / dense regime can be quantified by the hierarchy of scales α s α 2 s a 1/3 1, (1.1b) or equivalently, in term of the saturation momenta Q s,p and Q s,t of the projectile and target respectively, 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 qq pairs, 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äcker-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][43][44][45]. In this approach, observables such as forward double valence-quark production cross sections dσ qvqv [46], forward 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σ qvqvqvqv have 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 single-and double-pair production dσ qq and dσ (qq)(qq) 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 N c for the double-pair production cross section dσ (qq)(qq) 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 N c , 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 qq pairs 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σ (qq)G to 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 single-and 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 2 we 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 3 we compute the production cross section for a single qq pair 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 4 we proceed to calculate the new cross sections for the production of two sets of soft particles: double qq pair production in section 4.1, mixed qqG production in section 4.2, and double gluon production in section 4.3. The successful cross-check against the double-gluon production cross section JHEP02(2019)024 Figure 1. The light-front wave functions to radiate a soft qq pair 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 5 we 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 6 by reiterating the primary new results and exploring the many opportunities for phenomenological applications and further theoretical development which this work provides. In appendix A we provide details of the Gaussian color averaging used for the (semi-)dilute projectile, in appendix B we formulate some useful algebraic properties of Wilson lines in momentum space, and in appendix C we 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 with magnitudes v T ≡ |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.

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 qq pair as ψ 1 , ψ 2 , ψ 3 corresponding to the various time orderings of the scattering in the target fields. The term ψ 1 corresponds to scattering after the pair is created, ψ 2 to scattering after the gluon is emitted but before the pair is created, and ψ 3 to 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 where α ≡ k + k + +k + is the fraction of the pair longitudinal momentum carried by the quark, q is the center-of-mass transverse momentum of the qq pair (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 g from the definition of the wave functions. This corresponds to absorbing this coupling constant into the scale µ 2 defined 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]) where, as labeled in figure 1, x, y, and b are the final-state positions of the quark, antiquark, and valence quark, respectively, and u ≡ αx + (1 − α)y is the center-of-mass position of the qq pair (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, where we work in the A + = 0 light cone gauge. With (2.2) written this way, the Wilson line V b associated 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 b which 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 (V b t a ) → ρ a (b). (For another discussion of the translation between discrete and continuous charge distributions, see e.g. [51].) With this change of JHEP02(2019)024 notation, we can Fourier transform the buildling block (2.2) into momentum space to obtain where the momenta of the final-state quark and antiquark are k andk, respectively. Note that the Fourier factor for the valence quark cancels because its position b is 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 with q = k +k for brevity. We can combine all three diagrams by redefining the dummy integration variables, obtaining the compact form where the differences among the three diagrams are all encoded in the combined wave function  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 t a V † and the wave function Ψ, which is a matrix in the 2 × 2 spin space of the pair.

Gluon production amplitude
In comparison with (2.2) for the production amplitude of a soft qq pair in coordinate space, the corresponding amplitude to emit a soft gluon is illustrated in figure 2 and is given by where φ is the light-front wave function to radiate a soft gluon from a valence quark projectile. Here σ v and σ v are 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: (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 JHEP02(2019)024 emitted, and we have used the fact that the wave functions for the two time orderings differ by a minus sign (similar to ψ 1 + ψ 2 + ψ 3 = 0 for the quark pair case (2.1)). As we did in section 2.1, we can rewrite the second time ordering so that the valence quark scattering looks the same as the first one, and we can convert to the continuous charge density to write Comparing the gluon amplitude (2.12) with the pair amplitude (2.7), we note that the gluon has a specified color a in the final state, resulting in the pair-like structure V t b V † being contracted with 2t a to form a trace. Squaring the gluon amplitude gives where we denote the averaging over color fields of the projectile and target by · · · proj and · · · tgt , respectively. By tr c we denote a trace over color indices and by tr D we denote a trace over the 2 × 2 spin states of the wave function which averages over spins in the initial state and sums over spins in the final state. If we use the Fierz identity over the repeated color index a, we can combine the two traces into one, with the N c -suppressed term in the Fierz identity vanishing exactly by the Wilson line identity (B.5). This gives which has the same Wilson line structure as we would obtain by squaring the pair amplitude (2.7). Thus we can, without loss of generality, replace the gluon amplitude (2.12) with the equivalent expression that has the same form as the pair amplitude (2.7).
This similar structure appears to suggest a possible mapping between the pair amplitude (2.7) and the gluon amplitude (2.15). Comparing the two expressions, we see that in the limit k =k = 1 2 q, the Wilson line structure of the pair amplitude (2.7) can be cast in the same form as the gluon amplitude (2.15):

JHEP02(2019)024
The change of variablesk = q − k followed by k → −k + 1 2 q makes the comparison with (2.15) explicit, Thus we see that by mapping the pair wave function and setting the final-state quark and antiquark to have equal momenta, k =k = 1 2 q, we can map the pair amplitude (2.7) onto the gluon amplitude (2.15): Thus any cross section calculated for the production of multiple qq pairs via (2.7) can be mapped onto the cross section to instead produce gluons via (2.18) and (2.19). The existence of this mapping allows us to efficiently compute multiparticle production at mid rapidity, by first computing the production of various qq pairs and then mapping them systematically back to gluons. Aside from (2.19), the only other modifications to the cross section will be a change in the prefactor of 1 2(2π) 3 to reflect the changed number of finalstate particles and the exclusion of quark entanglement in the pair which has been mapped to a gluon. We will use this strategy in section 3.2 to obtain the single-inclusive gluon cross section from the pair cross section and in sections 4.2 and 4.3 to obtain the qqG and GG cross sections from the double-pair cross section.

Cross section for single-pair production
The inclusive cross section to produce a single soft qq pair with quark (antiquark) transverse momentum k (k) and rapidity y (ȳ) is simply related to the square of the amplitude (2.7): where dy = dk + k + and dȳ = dk + k + . To help keep the notation compact let us introduce the following notation for the differentials: (3.2c) We will also often exclude the underlines for the many transverse vectors when it is clear from context that they refer to 2-vectors rather than 4-vectors. Then squaring (2.7) and JHEP02(2019)024 performing the averaging as in appendix A, we straightforwardly obtain as illustrated in figure 3. Using the Gaussian averaging of the projectile with the assumption of Locality from (A.3), we obtain where the second argument of µ 2 vanishes in this case because the eikonal Wilson lines preserve the plus momentum. This feature will in general be violated when producing multiple pairs. The color trace is straightforward to simplify using the Fierz identity, yielding the compact expression where the fundamental dipole and quadrupole operators are defined in (A.1). Comparing this expression with the corresponding ones (32) and (37) from ref. [54] in coordinate space, we see that eq. (3.5) is far more compact. This is largely because the condensed amplitude (2.7) in momentum space combines the two time orderings "1" and "2" into a single form, allowing us to perform a single calculation for all time orderings, rather than requiring a sum over all the distinct time orderings for the pair emission. The single-pair cross section is then immediately given by combing eqs. (3.5) and (3.1).

Cross section for single-gluon production
As a cross-check and to illustrate the "gluonic mapping" of a qq pair (2.19), let us map (3.5) onto the cross section for single-inclusive gluon production. Applying (2.19) to (3.5) gives where we changed the prefactor from (3.1) to 1 2(2π) 3 to reflect the fact that there is now only one particle tagged in the final state. Changing variables to q ( , ) ≡ k ( , ) +k ( , ) and δk ( , At first glance, the mapping (2.19) doesn't seem to have accomplished very much. Even more alarmingly, this expression contains a quadrupole contribution, whereas the explicit gluon production cross section has only double-dipoles. The resolution is that there are hidden cancellations among the Wilson lines in momentum space, as described in appendix B. These cancellations in momentum space occur after integration over the dummy variables δk and δk , which do not couple to the wave functions or to µ 2 . Consider first the problematic quadrupole term: where we've used (B.5) to integrate over δk in the first two arguments and δk in the last two. The quadrupole term therefore sets q = q = q, but we immediately see from (2.18) that Φ(q; q) = 0, (3.9) so the quadrupole term vanishes identically in the "gluonic limit," as it must. This cancellation in the wave function is just a reflection of the fact that the sum of time orderings must vanish by the definition of the T -matrix. For the double-dipole term of (3.7), the situation is more subtle, since the shared momenta δk , δk live in opposite traces and This quantity, after integration over δk , δk , is just the Fourier transform to momentum space of a squared dipole amplitude. (Note that this is NOT the same thing as the square of the momentum-space dipole amplitude!) Using these results, the single-inclusive gluon cross section (3.7) becomes which is again an incredibly compact expression in momentum space. To complete the cross-check, let us insert the Fourier transform back to coordinate space, obtaining (3.12) From (2.18), it's straightforward to see that the Fourier transform of the wave function is where we have dropped the primes on any remaining integration variables in a given term. The last step is to insert the explicit gluon-emission wave functions (2.11), obtaining

JHEP02(2019)024
which agrees perfectly with eq. (38) of [54] and, after employing (B.6), with eq. (8.18) of [13]. This confirmation validates the cross-check of the "gluonic limit" (2.19) of a qq pair, as well as the dictionary (A.4) between the langauge of continous color charge densities and discrete valence quarks, since is the number of nucleons in the projectile. In the same way, we can next calculate the double-pair cross section and then use the mapping (2.19) to convert it into the qqG cross section and then into the GG cross section.

Cross section for double-pair production
For double-pair production, we consider a final state with two quarks k 1 , k 2 and two antiquarksk 1 ,k 2 . As emphasized in ref. [36], if the produced (anti)quarks have the same flavor, then we need to take care to explicitly antisymmetrize the full amplitude A full under the interchange of identical quarks and antiquarks to satisfy Fermi-Dirac statistics: Moreover, since in our case the two pairs are radiated independently, the amplitude is automatically symmetric under the interchange of both pairs: which can be seen because of the sum over color sources a in (2.7). As such, it is sufficient to only antisymmetrize the amplitude under the exchange of the antiquarks: with the elementary single-pair amplitudes being given by (2.7). Written explicitly, the Note that, as claimed, the amplitude is automatically symmetric under the exchange of the pairs (k 1 ,k 1 ) ↔ (k 2 ,k 2 ) if one also relabels the dummy momenta (k 1 ,k 1 ) ↔ (k 2 ,k 2 ) and JHEP02(2019)024 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 where the notation +(k 1 ↔k 2 ) 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.

Case 1: no fermion entanglement
Squaring (4.4) for topologies with no fermion entanglement, as in figure 4, leads directly to

JHEP02(2019)024
There are now 3 possible "contractions" of the source colors, obtained in terms of (A.3): 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: With this shorthand, the Wilson line tensor entering (4.6) is and we can straightforwardly compute the various contraction of (4.7): complete result for topologies with no fermion entanglement: Figure 6. Double-pair production "Pac Man" topologies with fermion entanglement, as calculated in eq. (4.12).

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 6 we have directly In the same way as (4.7), we form the 3 contractions of the projectile sources: = δ ab δ cd µ 2 (k 1 +k 1 + k 2 +k 2 , k + 1 +k

JHEP02(2019)024
Using the same shorthand notation as (4.8), the Wilson line tensor entering (4.12) is (4.14) and we can compute the various color contractions in the same way: (4.15c) Combining these back into (4.12) yields the complete result for topologies with fermion entanglement:

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 (qq) (qq) cross section at this order. These expressions are one of the primary results of this paper.

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 qq pair k 2 ,k 2 → 1 2 q 2 in 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 qqG are now distinguishable. Thus we need only take the limit of the topologies in (4.11), immediately obtaining: where as before we have changed variables into q 2 ≡ k 2 +k 2 and δk 2 ≡ 1 2 (k 2 −k 2 ), and similarly for q 2 and δk 2 .
After integration over δk 2 , δk 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 ( , ) 2 appears 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 (q 2 − q 2 ) or δ 2 (q 2 − q 2 ), 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): ≡ D 4 , 4 (p 1 , q 2 − q 2 , q 2 − q 2 , p 2 ; p 3 p 4 ). The result is the complete expression for the qqG cross section in momentum space, with the corresponding operators illustrated in figure 7: To our knowledge, eq. (4.19) represents the first complete calculation of qqG production in the CGC framework, and this compact form in momentum space comprises an exact solution at this order, at finite N c . 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.

Double-gluon production
As a final cross-check of the preceding calculations, let us use the mapping (2.19) on the remaining qq pair 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 k 1 ,k 1 → 1 2 q 1 in (4.19), adjust the count of 1 2(2π) 3 in the prefactor, and change integration variables to q 1 ≡ k 1 +k 1 and δk 1 ≡ 1 2 (k 1 −k 1 ), and similarly for q 1 and δk 1 . This gives Of the four interaction terms remaining in the braces, the third one (quadrupole trace) vanishes after integration over δk 1 δk 1 due 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):

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 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 C we 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 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.

JHEP02(2019)024
For the first term, the dependences on coordinates with subscript 1 and subscript 2 completely factorize: so that we can just multiply this result by the same expression under the exchange 1 ↔ 2.
For the second and third terms which share a set of "crossed" wave functions, the expressions do not cleanly factorize, but the delta functions generated do possess the following symmetry: The result of these delta functions, after dropping all primes, is to sum over all possible ways to take the squared quadrupole |D 4 | 2 (x 1 , y 1 , x 2 , y 2 ) and map x i → b i or y i → b i , generating a minus sign for each such replacement. Note that some of these permutations will cancel pairs of Wilson lines to generate squared dipoles and, in the last term, the unit operator. The result is The last step is to insert the explicit form of the gluon emission wave function (2.10) and to convert the squared fundamental dipole and quadrupole into the adjoint dipole and quadrupole using (B.6). The final expression for the double-gluon cross section is 2)) or may include heavy quarkonia such as (DD) (J/ψ) (center panel, eq. (5.6)). The partonic (cc) G cross section could hadronize similarly, but include a light hadron h produced by the gluon, as in (J/ψ) h (right panel, eq. (5.7)).
, so that the projectile is infinitely Lorentz-contracted into a delta function pancake at x − = 0. Our calculation here relaxes that assumption and is therefore sensitive to the 3D structure of the projectile. The successful cross-check (4.27) against ref. [34] serves as a validation of the preceding calculations based on the double-pair expressions (4.11) and (4.16).

Hadronization in the heavy flavor sector
To connect with the experimental measurements of interest, the partonic cross sections calculated above need to be convoluted with appropriate nonperturbative descriptions of hadronization. While there is an abundance of hadronic observables which are of interest, one of the cleanest analogs between the partonic and hadronic cross sections is in the heavy flavor sector. While for light hadron production, many partonic channels can all contribute with comparable weights, heavy flavor production is dominated by the production of heavy quark pairs [58,59]. We can characterize heavy flavor production in either the open or hidden heavy flavor sectors; here for specificity we will consider the production of D andD mesons in the open sector or J/ψ mesons in the hidden sector as illustrated in figure 8, and we will generally follow the methods of augmenting CGC calculations with hadronization prescribed in ref. [55]. The various hadronic cross sections presented in this section are the final primary result of this paper.
For production of open D mesons, we need to convolute the partonic cross sections with fragmentation functions following the standard techniques [60]. Thus from the singlecc cross section expressed in eqs. (3.1) and (3.5) we can directly compute the DD cross section via the convolution

JHEP02(2019)024
where the final-state momenta of the D andD are given by p µ D = zk µ and p μ D =zk µ , respectively in the limit where the hadrons can be considered massless. In this limit, the rapidities of the hadrons are equal to the rapidities of the quarks: y D = y and yD =ȳ, and the factors of 1/z 2 , 1/z 2 arise from the change of variables d 2 p D = z 2 d 2 k and d 2 pD = z 2 d 2k . The fragmentation functions D D/c and DD /c can be taken from e.g. ref. [61]. The lower limits of integration in (5.1) come from the eikonal approximation for the emission of the soft gluon from the projectile: k + ,k + 1 a P + a , where 1 a P + a is the average momentum per nucleon of a projectile with a nucleons whose total momentum is P + a . In principle the cc production mechanism calculated in eq. (3.5) is only one partonic channel that can contribute to final-state DD production, and one should generalize (5.1) by summing over all such partonic channels. However, as argued above, in the heavy flavor sector we expect the cc channel to make the dominant contribution to DD production. It is also important to note that the assumption of independent fragmentation of the c andc quarks relies on the partons being well-separated in phase space (either in transverse momentum or in rapidity) such that nonperturbative interference effects in hadronization can be neglected. Studies of the boundaries of independent fragmentation have been performed in the context of distinguishing target vs. current fragmentation regions in electron-proton collisions [62].
In the same way, we can convolute the (cc) (cc) cross section expressed in eqs. (4.5), (4.11), and (4.16) with fragmentation functions to obtain the four-particle hadronic cross section for DD DD production (left panel of figure 8): Here the final-state hadron momenta are related to the partonic ones by p µ D1 = z 1 k µ 1 , p μ D1 =z 1k µ 1 , p µ D2 = z 2 k µ 2 , and p μ D2 =z 2k µ 2 . Again, implicit in (5.2) is an assumption that the hadrons are well-separated in phase space, such that there is no interference during hadronization which mixes the pairs. This can be accomplished, for instance, by a rapidity gap 3 between the pairs y 1 ,ȳ 1 y 2 ,ȳ 2 and a large relative momentum between each corresponding D andD meson: |p D1 −pD 1 | T , |p D2 −pD 2 | T Λ QCD . Such kinematics are, in fact, appropriate for exploring the long-range correlations which may exist among the produced hadrons.
Likewise, we can also convolute the (cc) G cross section given in (4.19) with fragmentation functions to describe the three-particle correlations between a D meson, aD meson,

JHEP02(2019)024
and a light hadron h. That expression is given by Again, the final-state hadronic momenta are related to the partonic ones by p D = z k, pD =zk, and p h = z h q, and the usual caveats for independent fragmentation apply. For the production of heavy quarkonia such as the J/ψ, we will follow the procedure of [55] and employ the Improved Color Evaporation Model (ICEM) [63] to treat the hadronization of a partonic cc pair into a final-state J/ψ. For the hadronization of a single pair, we obtain the hidden heavy flavor analogue of (5.1): with M = (k +k) 2 the invariant mass of the cc pair, ∆k the relative momentum between the cc pair in its rest frame, and φ the corresponding angle. The ratio M m J/ψ is a conversion factor allowing the momentum P of the final J/ψ to differ from the cc center-of-mass momentum, and the other factors correspond to the Jacobian from the change of variables. The single nonperturbative parameter F J/ψ describes the probability for the cc pair to hadronize to the J/ψ by randomly emitting soft particles.
Similarly, we can translate the perturbative (cc)(cc) cross section from eqs. (4.5), (4.11), and (4.16) into a production cross section for two J/ψ mesons, writing the analogous expression Now M 1 is the invariant mass of the pair (k 1 ,k 1 ) and M 2 is the invariant mass of the pair (k 2 ,k 2 ), with the relative momenta ∆k 1 , ∆k 2 and angles φ 1 , φ 2 , respectively. Again, this JHEP02(2019)024 expression assumes a sufficient kinematic separation in phase space such that we do not have to consider interferences during hadronization. Last, we can consider hadronic cross sections which involve both the formation of a quarkonium state via the ICEM and fragmentation. The partonic (cc) (cc) cross section can hadronize into a J/ψ as well as an open DD pair (center panel of figure 8), which we write as Similarly, we can write the cross section for the hadronization of the partonic (cc) G state into a J/ψ meson in association with a light hadron h (right panel of figure 8): The above hadronic cross sections represent all possible ways for the partonic production of a (cc) (cc) state to hadronize into open DD pairs or J/ψ quarkonia, as well as all possible ways for a (cc) G state to hadronize into DD or J/ψ along with a light hadron h. These hadronic cross sections open the door to study a wealth of multi-hadron final states in pA and heavy-light ion collisions. Before concluding, a few comments are in order about the models and assumptions used for hadronization above. In the inclusive hadronization of a charm quark into a D meson or a gluon into a light hadron, we have employed collinear fragmentation functions. The use of collinear nonperturbative functions for fragmentation stands in contrast to the description of the initial state, which includes the effects of intrinsic transverse momentum arising from the nonperturbative wave functions of the projectile and target. This asymmetric treatment of the intrinsic transverse momentum scales is justified by the asymmetric nature of the semi-dilute/dense regime. The intrinsic transverse momentum of the dense target is characterized by the target saturation scale Q s,t and is resummed to all orders; JHEP02(2019)024 the intrinsic transverse momentum of the pojectile is characterized by the projectile saturation scale Q s,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 particularcarry 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. 4 In 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 qq pair, 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 cc cross section as done above, using a series of projection operators to select out the quantum numbers of the cc state appropriate for a given hadronization channel. Of these projection operators, the color projections onto singlet and octet cc states 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 JHEP02(2019)024 to quarkonium production for future work. It is interesting to note, however, that inclusive J/ψ production at small x in 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 3 S [8] 1 production channel.

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 (qq) (qq) (eqs. (4.5), (4.11), and (4.16)) and for one quark/antiquark pair plus a gluon (qq) 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 N c . 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 qq pair and the production amplitude for a gluon, which we used to obtain the (qq) G cross section from the (qq) (qq) 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 5 we 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: (DD) (DD) -eq. (5.2); (DD) (J/ψ) -eq. (5.6); (J/ψ) (J/ψ) -eq. (5.5); (DD) heq. (5.3); and (J/ψ) h -eq. (5.7), where h is 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 heavy-and heavy-light ion collisions. Genuine multiparticle production computed within the CGC framework makes it possible to self-consistently compute higher cumulants such as v 2 {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-N c approximation, we anticipate that it will open up 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 D 8 for 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: (qq) (qq) (qq). 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:(qq) (qq) (qq); (qq) (qq) G; (qq) G G; and G G G.
Last, we note that the expressions for hadronization in the heavy flavor sector we explore in section 5 can be significantly improved and extended. The heavy flavor sector is convenient as a justification for the assumption that a given hadron (like a D meson) in the final state is dominated by fragmentation from a given parton (like a c quark). 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.

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 · · · proj , separately from the averaging of the Wilson lines over the color fields in the target, which we denote · · · tgt or simply · · · when there is no ambiguity. The collision geometry is characterized by an impact parameter B between 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 bŷ For quantities that have already been averaged we drop the hat over the operator, writing e.g. D 2 (x, y) = D 2 (x, y) . When computing similar traces with adjoint Wilson lines, we will denote the operator with the superscript "adj", writing e.g.

JHEP02(2019)024
And we will maintain the same notation whether invoking Wilson lines in coordinate space or momentum space, writing e.g.D 2 (p, q) = 1 Nc tr[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 ⊥ 2 of the collision, and the quantum evolution with the running of this scale is given by the JIMWLK evolution equations [10][11][12] (or the large-N c analogue, 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 ∆y ij < 1 αs for the rapidity interval ∆y ij between any two particles i, j tagged at mid-rapidity.
For the average · · · proj 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 ρ ρ to fully specify the result of the averaging procedure: 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: δ ab µ 2 (q 1 − q 2 , q + 1 − q + 2 ) Locality δ ab (2π) 2 δ 2 (q 1 − q 2 ) µ 2 (q 2 1T , q + 1 − q + 2 ) 2D Transl. Inv. δ ab (2π) 2 δ 2 (q 1 − q 2 ) µ 2 (q + 1 − q + 2 ) Both. (A.3b) For averages with both of the color fields in the amplitude ρ ρ , the momentum of the second source is reversed: q 2 → −q 2 , and for averages with both of the color fields in the complex conjugate amplitude ρ * ρ * , the momentum of the first source is reversed: q 1 → −q 1 . 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 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: Clearly if we integrate over the shared momentum p, we generate a delta function that sets x = y and cancels the Wilson lines. The remaining integral over d 2 y then generates a new delta function: 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 so we can conclude that the two conditions are in fact equivalent: 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 N c suppressed: In going to momentum space, that additive constant becomes instead proportional to delta functions: conversion factor ( 4παs g 2 ) 2 is just a trivial difference of convention: we take our µ 2 to include the gluon-emission coupling g 2 , while they insert the coupling explicitly. We do, however, differ by the numerically significant prefactor 16 [2(2π) 3 ] 2 ≈ 6.5 × 10 −5 which 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.