Parton model description of multiparticle azimuthal correlations in $pA$ collisions

TL;DR

This paper extends a parton model to analyze multiparticle azimuthal correlations in $pA$ collisions, challenging hydrodynamic interpretations and exploring the effects of multiple scatterings, saturation, and transverse momentum dependence.

Contribution

It introduces a detailed analysis of azimuthal correlations in the parton model, including $p_ot$ dependence, $N_c$ effects, and an algorithm for Wilson line correlators, expanding prior work and providing new insights.

Findings

Coherent multiple scattering is crucial for positive $v_2 ext{ }\{4\}$.

Azimuthal correlations depend on the number of color domains and $p_ot$ window.

Extended the analysis to symmetric cumulants and $N_c$ dependence.

Abstract

In arXiv:1705.00745, an initial state "parton model" of quarks scattering off a dense nuclear target was shown to qualitatively reproduce the systematics of multiparticle azimuthal anisotropy cumulants measured in proton/deuteron-nucleus ($pA$) collisions at RHIC and the LHC. The systematics included i) the behavior of the four-particle cumulant $c_2\{4\}$, which generates a real four-particle second Fourier harmonic $v_2\{4\}$, ii) the ordering $v_2\{2\}>v_2\{4\}\approx v_2\{6\}\approx v_2\{8\}$ for two-, four-, six-, and eight-particle Fourier harmonics, iii) the behavior of so-called symmetric cumulants $\text{SC}(2,3)$ and $\text{SC}(2,4)$. These features of azimuthal multiparticle cumulants were previously interpreted as a signature of hydrodynamic flow; our results challenge this interpretation. We expand here upon our previous study and present further details and novel results on the saturation scale and transverse momentum ($p_\perp$) dependence of multiparticle azimuthal correlations. We find that the dependence of $v_2\{2\}$ and $v_2\{4\}$ on the number of color domains in the target varies with the $p_\perp$ window explored. We extend our prior discussion of symmetric cumulants and compute as yet unmeasured symmetric cumulants. We investigate the $N_c$ dependence of $v_2\{2\}$ and $v_2\{4\}$. We contrast our results, which include multiple scatterings of each quark off the target, to the Glasma graph approximation, where each quark suffers at most two gluon exchanges with the target. We find that coherent multiple scattering is essential to obtain a positive definite $v_2\{4\}$. We provide an algorithm to compute expectation values of arbitrary products of the "dipole" lightlike Wilson line correlators.

Equations

$${⟨q(\mathbf{q}{)}_{out}|q(\mathbf{p}{)}_{in}⟩}_{A^{-}}=\int d^2{\mathbf{x}}_{\perp }\left[U(x_{\perp })-1\right]e^{i(\mathbf{q}-\mathbf{p})\cdot {\mathbf{x}}_{\perp }}\sim \mathcal{M}(p,q),$$

$$U({\mathbf{x}}_{\perp })$$

$$⟨\frac{dN}{d^2\mathbf{p}}⟩$$

$$W_q(\mathbf{b},\mathbf{k})=\frac{1}{{\pi }^2}{e}^{-|\mathbf{b}{|}^2/B_p}{e}^{-|\mathbf{k}{|}^2{B}_p},$$

$$D(x,y)=\frac{1}{N_c}\text{Tr}\left[U(x)U^{†}(y)\right]$$

$$⟨\frac{dN}{d^2\mathbf{p}}⟩$$

$$⟨\frac{d^{m}N}{d^2{\mathbf{p}}_1\cdots d^2{\mathbf{p}}_{\mathbf{m}}}⟩\equiv ⟨\frac{dN}{d^2{\mathbf{p}}_1}\cdots \frac{dN}{d^2{\mathbf{p}}_{\mathbf{m}}}⟩,$$

$$W({\mathbf{b}}_1,{\mathbf{k}}_1,...,{\mathbf{b}}_{\mathbf{m}},{\mathbf{k}}_{\mathbf{m}})=W_q({\mathbf{b}}_1,{\mathbf{k}}_1)\cdots \cdots W_q({\mathbf{b}}_{\mathbf{m}},{\mathbf{k}}_{\mathbf{m}}).$$

$$⟨\frac{d^{m}N}{d^2{\mathbf{p}}_1\cdots d^2{\mathbf{p}}_{\mathbf{m}}}⟩$$

$$[D_{\mu },F^{\mu \nu }]={\delta }^{\nu -}\rho (\mathbf{x},x^{+}),$$

$$⟨{\rho }_a({\mathbf{x}}_{\perp },x^{+}){\rho }_b({\mathbf{y}}_{\perp },y^{+})⟩={\delta }_{ab}\delta (x^{+}-y^{+}){\delta }^{(2)}({\mathbf{x}}_{\perp }-{\mathbf{y}}_{\perp }){\mu }^2(x^{+}),$$

$$A_a^{-}({\mathbf{x}}_{\perp },x^{+})=g\int d^2{\mathbf{z}}_{\perp }G({\mathbf{x}}_{\perp }-{\mathbf{z}}_{\perp }){\rho }_a({\mathbf{z}}_{\perp },x^{+}),G({\mathbf{x}}_{\perp })=\int \frac{d^2{\mathbf{k}}_{\perp }}{(2\pi {)}^2}\frac{e^{i{\mathbf{k}}_{\perp }\cdot {\mathbf{x}}_{\perp }}}{|{\mathbf{k}}_{\perp }{|}^2},$$

$$g^2⟨A_a^{-}({\mathbf{x}}_{\perp },x^{+})A_b^{-}({\mathbf{y}}_{\perp },y^{+})⟩$$

$$L({\mathbf{x}}_{\perp },{\mathbf{y}}_{\perp })$$

$$⟨D\cdots D⟩=\mathcal{T}\mathcal{N}.$$

$$t_{ij}^a{t}_{kl}^a=\frac{1}{2}{\delta }_{il}{\delta }_{jk}-\frac{1}{2N_c}{\delta }_{ij}{\delta }_{kl},$$

$$\displaystyle\mathcal{T}=\text{exp}\Big{(}{-\frac{C_{F}}{2}\displaystyle\sum_{i=1}^{4}\big{(}L_{x_{i}x_{i}}+L_{\bar{x}_{i}\bar{x}_{i}}\big{)}}\Big{)},$$

$$U({\mathbf{x}}_{\perp })$$

$$⟨D({\mathbf{x}}_{\perp },{\bar{\mathbf{x}}}_{\perp }{)}_U⟩$$

$$⟨D_{x_1{\bar{x}}_1}{D}_{x_2{\bar{x}}_2}{⟩}_U$$

$${\left(\begin{array}{c}⟨D_{x_1{\bar{x}}_1}{D}_{x_2{\bar{x}}_2}⟩\\ ⟨Q_{x_1{\bar{x}}_2{x}_2{\bar{x}}_1}⟩\end{array}\right)}_U={\left(\begin{array}{cc}{\alpha }_{x_1{\bar{x}}_1{x}_2{\bar{x}}_2}& {\beta }_{x_1{\bar{x}}_2{x}_2{\bar{x}}_1}\\ {\beta }_{x_1{\bar{x}}_1{x}_2{\bar{x}}_2}& {\alpha }_{x_1{\bar{x}}_2{x}_2{\bar{x}}_1}\end{array}\right)}_U{\left(\begin{array}{c}⟨D_{x_1{\bar{x}}_1}{D}_{x_2{\bar{x}}_2}⟩\\ ⟨Q_{x_1{\bar{x}}_2{x}_2{\bar{x}}_1}⟩\end{array}\right)}_V$$

$$⟨D_{x_1{\bar{x}}_1}{D}_{x_2{\bar{x}}_2}{D}_{x_3{\bar{x}}_3}{D}_{x_4{\bar{x}}_4}{⟩}_U$$

$${\vec{\mathcal{N}}}_n={\text{a}}_n{\vec{\mathcal{N}}}^{(\text{a})}+{\text{b}}_n{\vec{\mathcal{N}}}^{(\text{b})}+...+{\text{x}}_n{\vec{\mathcal{N}}}^{(\text{x})},$$

$${\vec{\mathcal{N}}}_{n+1}=\mathbf{M}{\vec{\mathcal{N}}}_n.$$

$$\left(\begin{array}{c}{\text{a}}_{n+1}\\ {\text{b}}_{n+1}\\ ...\\ {\text{x}}_{n+1}\end{array}\right)=\mathbf{M}\left(\begin{array}{c}{\text{a}}_n\\ {\text{b}}_n\\ ...\\ {\text{x}}_n\end{array}\right).$$

$$\vec{\mathcal{N}}=\sum _{n=0}^{\infty }\frac{1}{n!}{\vec{\mathcal{N}}}_n=\sum _{n=0}^{\infty }\frac{1}{n!}{\mathbf{M}}^n{\vec{\mathcal{N}}}_0=e^{\mathbf{M}}{\vec{\mathcal{N}}}_0.$$

$$\mathcal{N}=(1,1/N_c,..,1/N_c,1/N_c^2,...,1/N_c^2,1/N_c^3,...,1/N_c^3)\vec{\mathcal{N}}.$$

$$⟨DDDD⟩$$

$$U({\mathbf{x}}_{\perp })=\mathcal{P}\exp \left(-ie\int d{z}^{+}{A}^{-}({\mathbf{x}}_{\perp },z^{+})\right).$$

Full text

I Introduction

In Ref. Dusling:2017dqg , we demonstrated that a simple initial state parton model gives rise to many of the features of multiparticle azimuthal correlations observed experimentally in small collision systems at both RHIC at BNL and the LHC at CERN. In this model, collinear quarks from the projectile scatter coherently off color sources of the size of the inverse of the saturation scale $Q_{s}$ in the nuclear target; the scattered quarks are found to be collimated in their relative azimuthal angles. We observed crucially that one obtains a negative value for $c_{2}\{4\}$, the four-particle azimuthal anisotropy cumulant. This results in a positive definite four-particle Fourier coefficient $v_{2}\{4\}$. We demonstrated further, in a simpler Abelian version of our model, that one obtains the ordering of $m$-particle second Fourier harmonics $v_{2}\{2\}>v_{2}\{4\}\approx v_{2}\{6\}\approx v_{2}\{8\}$. Both of these features were previously believed to be unique signatures of collectivity arising from the hydrodynamic flow of quark-gluon matter. Not least, we demonstrated that so-called symmetric cumulants (mixed four-particle cumulants of different Fourier harmonics) computed in this parton model display the same qualitative features as experimental measurements of symmetric cumulants. We note that symmetric cumulants were designed to probe correlations and fluctuations arising from the hydrodynamic response to different harmonics of the azimuthal structure of the initial geometry.

While the hydrodynamic description of the flow of quark-gluon matter is likely valid in the larger, more central, collisions of heavy-ions ($AA$), the applicability of this description to more peripheral $AA$ collisions, and to $pA$ and proton-proton ($pp$) collisions, is less clear Bzdak:2013zma ; Niemi:2014wta . Our results in Ref. Dusling:2017dqg are therefore a strong hint that the stated measures of hydrodynamic collectivity are not robust in their own right without further corroboration from other distinct measures of collectivity. An example of the latter is the strong jet quenching that is seen in central heavy-ion collisions at RHIC and the LHC Adams:2003im ; Adams:2003kv ; Abelev:2012hxa ; CMS:2012aa . In contrast, jet quenching is either small or absent in peripheral $AA$ collisions and in $pA$ collisions Adams:2003im ; Adam:2016dau ; McGlinchey:2017esf .

We will elaborate in this paper on the parton model description introduced in Ref. Dusling:2017dqg , which extends previous work on two-particle azimuthal correlations discussed in Lappi:2015vha ; Lappi:2015vta . A novel feature is the development of a general algorithm (based on the framework in Blaizot:2004wv ) to compute expectation values of multi-dipole correlators. These objects encode the physics of multiple eikonal scattering of quarks on a colored target. In particular, the dipole operator is the trace over the product of a lightlike Wilson line appearing in the quark production amplitude at a given transverse position, with its conjugate transpose appearing in the complex-conjugate amplitude at a different spatial location, normalized by the number of colors $N_{c}$. We will present a systematic study of azimuthal cumulants and Fourier harmonics as a function of the target saturation scale $Q_{s}$ and the transverse momentum. Additionally, we will perform a systematic study of the $N_{c}$ dependence of observables. In particular, we will point to key similarities and differences between the non-Abelian and Abelian versions of the model. We will also make predictions of yet to be measured symmetric cumulants for higher order Fourier harmonics.

Expectation values over the dipole correlators are computed in the McLerran-Venugopalan (MV) model McLerran:1993ka ; McLerran:1993ni . This model includes coherent multiple scattering of the quarks in the projectile off the nuclear target. If we include at most two scatterings of the quarks, corresponding to the expansion of the Wilson lines to lowest nontrivial order, the expectation values correspond to the Glasma graph approximation. This approximation is applicable for $p_{\perp}>Q_{s}$. A model including quantum evolution of the Glasma graphs in the Color Glass Condensate (CGC) framework Iancu:2003xm ; Gelis:2010nm was previously applied to successfully describe key features of azimuthal correlations for $p_{\perp}\geq Q_{s}$ Dumitru:2008wn ; Dumitru:2010iy ; Dusling:2012iga ; Dusling:2012cg ; Dusling:2012wy ; Dusling:2013qoz ; Dusling:2015rja . We show that the Glasma graph correlators only produce positive values of the four-particle cumulant $c_{2}\{4\}$ and therefore do not correspond to a real $v_{2}\{4\}$. This result demonstrates that coherent multiple scattering, which is significant for $p_{\perp}\leq Q_{s}$, is an essential ingredient for the collectivity seen in our initial state framework.

The organization of the paper is as follow. We begin with the setup of our model in Sec. II. In Sec. III, we discuss the algorithm for the computation of multi-dipole correlators. Results of our computations are presented in Sec. IV. In Sec. V, we discuss the dependence of these results on the relative separations of the quarks in the projectile, on the number of color domains for varying $p_{\perp}$ windows and on the number of colors $N_{c}$. We contrast our results to those in the Glasma graph approximation. We briefly discuss the rapidity dependence on correlations in our model. In Sec. VI, we conclude and discuss possible future directions of research. Details of the Glasma graph computations are presented in the Appendix.

II Eikonal quark scattering from a nuclear target

We will discuss in this section a simple parton model description of proton-nucleus collisions. The incoming projectile consists of a collection of independent, nearly collinear, quarks that scatter off a dense nuclear target. Spatial correlations within the classical field of the nucleus imprint themselves on the quarks as they scatter, resulting in nontrivial momentum space correlations between the originally uncorrelated quarks. These include correlations in their relative azimuthal angles.

We begin by considering the scattering of a fermion off a classical background field in the high energy limit Bjorken:1970ah ; Mueller:1994gb . The forward scattering amplitude for a fixed background field $A^{-}$ can be expressed as Dumitru:2002qt

[TABLE]

where

[TABLE]

is the Wilson line in the fundamental representation at a transverse position $\mathbf{x}_{\perp}$ and $\mathcal{P}$ denotes path ordering in the lightcone variable $x^{+}$. The $-1$ in Eq. (1) removes the “no scattering” contribution wherein a quark passes through the target nucleus without having its color rotated by an Eikonal phase. As the incoming partons all have transverse momentum of order $\Lambda_{\rm QCD}$, and we are interested in $|\mathbf{p}|\gg\Lambda_{\rm QCD}$, we will ignore this contribution in what follows.

The transverse spatial distribution of collinear quarks with transverse momenta $\mathbf{k}$ in the projectile is represented by the Wigner function $W_{q}(\mathbf{b},\mathbf{k})$ Lappi:2015vha ; Lappi:2015vta . The single inclusive distribution within this model can be expressed as

[TABLE]

where the expectation value denotes an average over fields $A^{-}$ in the target, as for instance given by the MV model. For simplicity, we assume the Wigner function has the Gaussian form

[TABLE]

where both the transverse momentum and spatial location of the quarks is determined by a single nonperturbative scale $B_{p}$. Unless otherwise mentioned, we will fix $B_{p}=4~{}\text{GeV}^{-2}$ Kowalski:2006hc , obtained from dipole model fits to HERA deep inelastic scattering data. We will discuss later how our results are affected by variations in the value of $B_{p}$.

In Eq. (6), the function

[TABLE]

denotes the dipole operator. This operator encodes all orders in multiple gluon exchanges, as we will explicitly see in the calculations in Sec. III. Performing the integration over the incoming quark momenta $\mathbf{k}$, Eq. (6) can be simplified to read,

[TABLE]

The above framework can be extended to multiparticle production. For $m$ incoming quarks in the projectile the $m$-particle inclusive spectrum can be expressed as

[TABLE]

where the expectation value denotes an average over classical configurations of the target in a single event and over all events. Since each of the single-particle distributions inside the average here is a gauge-dependent functional of the classical field, we caution the reader that these distributions are qualitatively different from the gauge invariant single-particle distributions employed in hydrodynamic computations. No such simple product of gauge invariant distributions can be written in our case; indeed, as discussed at length in the Appendix, the Feynman diagrams corresponding to Eq. (7) are quantum interference diagrams.

On the other hand, the quarks comprising the projectile are uncorrelated, with the $m$-quark Wigner function of the projectile factorizing into a product of single-quark Wigner functions,

[TABLE]

Diagrammatically, this can be represented, as shown in Fig. (1), as the multiple scattering between different quarks in the amplitude and complex conjugate amplitude and the target nucleus. In the strict dilute-dense limit in which we work, these Wigner functions are gauge invariant distributions of which the product form is assumed to survive color averaging. From Eq. (7) and Eq. (8), we arrive at the following compact form for the $m$-particle inclusive spectra:

[TABLE]

Even though the above expression has a factorized form, it is highly nontrivial. Multiparticle correlations are generated via the expectation value over the classical fields of the target. A primary focus of this work will be computing the correlation between four particles. In this case, the expectation value is over a product of four-dipole operators, each of which, as noted, is a trace of two lightlike Wilson lines. The resulting expectation value is a function of eight transverse coordinates: four coordinates in the amplitude and four in the complex-conjugate amplitude. This expectation value will be evaluated without approximation in the MV model. We will see that large $N_{c}$ approximations and perturbative expansions (such as the Glasma graph approximation) to such correlators are insufficient to capture the systematics of $pA$ data.

A shortcoming of our model is the oversimplified nature of the projectile. At high energies, gluon radiation dominates the small-$x$ component of the proton’s wavefunction. These high parton densities become apparent when $Q_{s,T}/p_{T}\gtrsim 1$; saturation model fits to HERA data conservatively suggest that these effects become non-negligible around $x=0.01$ Rezaeian:2012ji . However, depending on the transverse momentum range studied, the qualitative features we observe could persist to smaller values of $x$. Furthermore, as the rapidity separation between quarks becomes larger than $\Delta y\gtrsim 1/\alpha_{S}$, quantum corrections will result in a decorrelation between partons. A more quantitative discussion of the rapidity dependence is discussed in Sec. V.5.

Quantum evolution will clearly break the factorized form of the Wigner function used in Eq. (8). Furthermore, since gluons would dominate the scattering process, their interactions with the target would be represented as adjoint Wilson lines. Multiparticle production, in this dense-dense limit, has been addressed in previous work Gelis:2008ad ; Dusling:2009ni ; Lappi:2009xa ; Schenke:2015aqa ; McLerran:2016snu ; Kovner:2016jfp . Multiparticle distributions can be obtained by solving the classical Yang-Mills equations in the presence of lightlike color sources corresponding to the projectile and the target. These source densities are each drawn from functional distributions of color charges, the evolution with energy of which is described by the JIMWLK equations JalilianMarian:1997gr ; Iancu:2000hn . However, due to the numerical complexity of the simulations, this has been restricted thus far to two-particle correlations Lappi:2009xa ; Schenke:2015aqa . The success of our simple model in explaining many of the collective signatures seen in light–heavy-ion collisions should stimulate further development of classical Yang-Mills simulations.

One may note that Eq. (9) has the structure of an expectation value of a product of functions. If one interpreted these functions as “single-particle” distributions the form of Eq. (9) would be markedly similar to a hydrodynamic framework Gronqvist:2016hym . One may then conjecture that the results we show for $v_{2}\{4\}$ are simply a consequence of the functional form of Eq. (9). This turns out not to be the case. In our discussion of coherent multiple scattering versus Glasma graphs, we will observe that, while both can be expressed in single-particle product form, one obtains negative four-particle azimuthal cumulants in the former case and positive valued cumulants in the latter111We thank Jean-Yves Ollitrault for a useful discussion on this point..

III Expectation values of multi-dipole correlators in the MV model

In this section, we will compute the expectation value of four dipole operators in the MV model McLerran:1993ka ; McLerran:1993ni of a nucleus at high energies. Although we will only make use of dipole operators, the results also generate all allowed expectation values of eight Wilson lines at no additional computational cost. The algorithm presented here to compute expectation values of lightlike Wilson line correlators can in principle be extended to higher-point functions.

In the MV model, classical gauge fields are described by solutions of the classical Yang-Mills equations,

[TABLE]

where $\rho$ denotes the classical color charge density in the nucleus. It is determined from the random Gaussian distribution satisfying

[TABLE]

where $\mu^{2}$ is the squared color charge density per unit area. The above two-point function can also be recast in terms of the gauge fields $A^{-}_{a}$ using

[TABLE]

where $G(\mathbf{x}_{\perp})$ is the free gluon propagator in two dimensions Gelis:2008rw . One then obtains,

[TABLE]

The integral over the two-dimensional propagator is formally divergent and must be regulated at the nonperturbative scale $\Lambda_{\rm QCD}$ with the result

[TABLE]

Expectation values of multiple dipoles can be computed by expanding the path ordered exponential within the Wilson line to second order in the gauge field. All possible pairwise contractions of the gauge fields are evaluated using Eq. (13). The result can be re-exponentied resulting in an expression valid to all orders in the gauge field Kovner:2001vi ; Fujii:2002vh . There is also a formally equivalent graphical method Blaizot:2004wv ; Dominguez:2008aa , which we will refer to when helpful.

When expanding out the path ordered exponential in the Wilson line, each gauge field represents a single gluon exchange with the target. Each $\langle AA\rangle$ contraction is therefore equivalent to considering two gluon exchanges between two quarks (represented as Wilson lines) and the target nucleus. Gluon exchanges can occur twice on the same quark, or on different quarks, and can also act on anti-quarks in the conjugate amplitude.

Following Blaizot:2004wv , we will denote as $\mathcal{T}$ contributions arising from two gluon exchanges with the same quark and $\mathcal{N}$ as two-gluon exchange among different quarks. Since exchanges between the same Wilson line in $\mathcal{T}$ are color singlets, these contributions can be considered separately from the $\mathcal{N}$ contributions. The final result for the expectation value of $m$-dipole operators can be expressed as a product of the two contributions,

[TABLE]

III.1 Tadpole contribution

Gluon exchanges between the same Wilson line comprise the tadpole contribution $\mathcal{T}$. We take as our starting point eight Wilson lines having transverse positions $\mathbf{x}_{1},\mathbf{\bar{x}}_{1},\dots,\mathbf{x}_{4},\mathbf{\bar{x}}_{4}$, where the $\mathbf{x}$ positions refer to quarks and the $\mathbf{\bar{x}}$ positions refer to anti-quarks. These dipoles must be connected in such a way to preserve the flow of color. Namely, quarks can only be connected to anti-quarks and vice versa. Without loss in generality, we connect at $x^{+}=-\infty$ quarks at $\mathbf{x}_{i}$ with anti-quarks at $\mathbf{\bar{x}}_{i}$, as shown in Fig. (2).

This allows us to unambiguously define what we mean as a four-dipole configuration for the given positions222Physically however these positions are completely arbitrary., as shown in Fig. 4(a). However, we can consider the $x^{+}=+\infty$ ends to be initially open ended. To evaluate the tadpole contribution, we draw a gluon connecting at two arbitrary $x^{+}$ points on the same Wilson line of Fig. (2). We then invoke the Fierz identity,

[TABLE]

which is given graphically by Fig. 3.

Performing this Fierz-ing, we see that the first term gives us a closed loop, which just is the trace of the identity matrix. This gives a factor of $N_{c}$ and the dipole again. The second term is just the original dipole. The result is then the original dipole multiplied with the Casimir factor $C_{F}=\frac{N_{c}^{2}-1}{2N_{c}}$.

For the rest of the result, we must calculate pair-wise contractions of the gauge field $A^{-}_{a}$. As noted, for the MV model, these satisfy Eq. (13). Every gluon exchange between the same (anti)-quark (same transverse position Wilson line) results in a factor of $-\frac{C_{F}}{2}L_{x_{i}x_{i}}$, where $L$ was defined in Eq. (14), and the $1/2$ is on account of the fact that the two ends of the correlator are ordered in $x^{+}$ because they belong to the same Wilson line. The negative sign is from connecting a (anti)-quark with a (anti)-quark.

Summing $n$ such exchanges to the Wilson lines, the tadpole contribution can be expressed as

[TABLE]

where we now denote the transverse position arguments as subscripts for readability and differentiate between quarks in the amplitude ($x_{i}$) and anti-quarks in the complex-conjugate amplitude with an over-bar ($\bar{x}_{i}$). This tadpole term is a color singlet and commutes with the terms we will derive next.

III.2 Gluon exchange contribution

We shall now consider gluon exchanges between different Wilson lines. Our starting point is again the configuration shown in Fig. (2). Closing the ends of the dipoles while preserving the color flow, we find that there are five distinct topologies. Using Fierz ordering as we did previously for the tadpole contribution, a gluon exchange between different Wilson lines in one topology can transform it into a different topology.

We now show how one obtains the five distinct topologies shown in Fig. 4. As a concrete example, we start with the four closed dipoles in Fig. 4(a). Consider a gluon exchange between $\bar{x}_{3}$ and $x_{4}$. From the two terms that result from Fierz ordering, as shown in Fig. 3, the first term gives two dipoles and a quadrupole, as depicted in Fig. 4(b), with a factor of $\frac{1}{2}$. A quadrupole, as depicted, is a distinct topological configuration corresponding to the trace over the product of Wilson lines at two distinct transverse spatial positions in the amplitude and at two such positions in the complex-conjugate amplitude. The second term from Fierz-ing just returns the original four-dipole configuration shown in Fig. 4(a), but now with the Fierz factor $\frac{1}{2N_{c}}$. Now taking this dipole-dipole-quadrupole configuration, consider further an exchange between $\bar{x}_{1}$ and $x_{2}$. This creates a quadrupole-quadrupole topology from the first term in the Fierz-ing, as depicted in Fig. 4(c), and likewise the structure Fig. 4(b) from the second Fierz term, both terms with the appropriate Fierz prefactors.

If instead we had started with the dipole-dipole-quadrupole configuration and considered a gluon exchange between $\bar{x}_{2}$ and $x_{3}$, the first Fierz term would have resulted in a dipole-sextupole configuration, as depicted in Fig. 4(d). A sextupole, as depicted in Fig. 4(d), corresponds to a trace over the product of Wilson lines at three spatial positions in the amplitude and three in the complex-conjugate amplitude. An exchange between $\bar{x}_{1}$ and $x_{2}$ in this dipole-sextupole topology results in an octupole (trace of eight Wilson lines–four in the amplitude and the other four in the complex-conjugate amplitude), depicted in Fig. 4(e), for the first Fierz term. The second term, as for the previous cases, gives back the original configuration, the dipole-sextupole one, with the appropriate Fierz prefactor. Thus, we see that multiple gluon exchanges continually generate, with each additional exchange, transitions between five topologically distinct configurations: (a) four-dipole, (b) dipole-dipole-quadrupole, (c) quadrupole-quadrupole, (d) dipole-sextupole, (e) octupole.

There is a transverse coordinate permutation degeneracy to these diagrams as well. The four-dipole is the only topology without this permutation degeneracy. For the dipole-dipole-quadrupole topology, there are $\binom{4}{2}=6$ possible permutations to close the $x^{+}=+\infty$ side of the eight pairwise connected Wilson lines we introduced previously. Similarly, there are $\frac{1}{2}\binom{4}{2}=3$ quadrupole-quadrupole permutations, $2\binom{4}{3}=8$ dipole-sextupole permutations, and $3!=6$ octopole permutations. As a sanity check, since we have considered eight Wilson lines connected pairwise at one end ($x^{+}=-\infty$), the sum of these permutations degeneracies agrees with the total $4!=24$ possible different contractions of the $x^{+}=+\infty$ side, which is dictated from the fact that quark (anti-quark) Wilson lines can only connect with anti-quark (quark) Wilson lines.

It should be clear from our discussion that we have a closed system of 24 configurations whereby each of these are transformed, through gluon exchanges and Fierz-ing, into other configurations in this system. We can express the 24 possible configurations as elements of a basis characterizing the four-dipole system. Starting with this initial condition, where all other configurations are set to zero, one can construct a 24 by 24 transformation matrix $\mathbf{M}$ that transforms one set of basis elements to another with each gluon exchange and subsequent Fierz-ing. We can then deduce, either analytically Kovner:2001vi ; Fujii:2002vh ; Fukushima:2007dy or diagramatically Blaizot:2004wv ; Dominguez:2008aa , what factors (via Fierz) are picked up in going from a basis element $\alpha$ to basis element $\beta$, which then define the elements $M_{\alpha\beta}$ of the matrix.

To understand how one fills in the arrays of this matrix, consider a path ordered exponential for the Wilson line

[TABLE]

where we expanded out the last infinitesimal slice in rapidity, and $V(\mathbf{x}_{\perp})$ is a redefinition of the original Wilson line excluding this last infinitesimal slice. Substituting this last expression in the dipole operator, we obtain

[TABLE]

Here we have made use of the locality in rapidity of correlators in the MV model.

Using Eq. (13), we can express $\langle D_{U}\rangle$ in the l.h.s in terms of $\langle D_{V}\rangle$ alone on the r.h.s.. Iterating this expression for each slice in rapidity, one obtains the exponentiated expression

[TABLE]

While this is a simple example, an identical procedure can be followed for multiple dipoles, or any higher order configuration, resulting from an initial four-dipole configuration. Another example is the well-known two-dipole result Blaizot:2004wv ; Dominguez:2008aa ; Kovner:2001vi ; Fujii:2002vh . In this case, following the procedure outlined above, it is straightforward to show that

[TABLE]

where $Q={\rm Tr}\left(V(x_{1})V^{\dagger}(\bar{x}_{2})V(x_{2})V^{\dagger}(\bar{x}_{1})\right)/N_{c}$ is the quadrupole configuration. Since the only two configurations for two dipoles are the dipole-dipole and quadrupole, we can then write this as a matrix equation:

[TABLE]

where $\alpha$ and $\beta$ are simple functions of $L(\mathbf{x}_{\perp},\mathbf{y}_{\perp})$ given in Ref. Blaizot:2004wv .

For a further example, if we consider four dipoles, then it is only possible via one gluon exchange to get to the six possible dipole-dipole-quadrupole configurations or stay in the four-dipole configuration. One obtains

[TABLE]

We can repeat this for all topologies to obtain a similar matrix as for the two-dipole case.

For a small number of dipoles, this procedure can be carried out efficiently by hand, and the eigenvalues of the matrix can be even computed analytically. However for larger numbers of dipoles, this becomes cumbersome. For example, for the four dipoles we have been considering, we must compute $7!$ gluon exchanges for each of the 24 basis elements, totaling 120,960 computations. Fortunately, since the algorithm suggested by our exercises is quite straightforward, it is very tractable to determine the elements of the 24 by 24 matrix and compute their eigenvalues on a computer. This work made extensive use of the GNU Scientific Library Gough:2009:GSL:1538674 and the EXPOKIT software package EXPOKIT . The matrix is illustrated schematically by Fig. 5; as it indicates, the elements of the matrix are relatively sparse but must nevertheless be diagonalized numerically. Our algorithm can be generalized to a larger number of dipoles of more complex topologies. It is therefore potentially useful for computations of many-body final states in high energy QCD where $n$-tupoles of lightlike Wilson lines are ubiquitous.

As a final remark on this matrix computation, we note that in the large $N_{c}$ limit this problem becomes more tractable Dominguez:2012ad ; Shi:2017gcq . In Ref. Dominguez:2012ad , it was shown that for large $N_{c}$ only dipoles and quadrupoles contribute to high energy QCD processes. However the azimuthal cumulants themselves vanish at large $N_{c}$, so in order to compute these quantities a finite $N_{c}$ calculation is necessary.

III.3 Result for product of four-dipole correlators

Now that we have explained how to calculate the matrix for one gluon exchange at a single slice in rapidity, we can compute how the basis vector space of configurations evolves after an infinite number of gluon exchanges. More generally, after $n$ gluon exchanges, we have the vector

[TABLE]

where $\vec{\mathcal{N}}^{(a)}$ through $\vec{\mathcal{N}}^{(x)}$ refer to individual basis vectors in the 24-dimensional basis space. This then evolves according to

[TABLE]

One can rewrite this expression as the matrix equation

[TABLE]

The setup of the computation thus far is general. However, as previously stated, we wish to begin with an initial condition that is the closed four-dipole configuration. This is done by setting the initial condition $\vec{\mathcal{N}}_{0}=\vec{\mathcal{N}}^{(a)}$ (the four-dipole configuration). Since we know how to compute additional gluon exchanges from Eq. (26), we need to multiply our initial condition $\vec{\mathcal{N}}_{0}=\vec{\mathcal{N}}^{(a)}$ by $\mathbf{M}$ $n$ times. This is, however, just a compact way to write all possible configurations with appropriate factors after $n$ gluon exchanges from the starting point of four closed dipoles. The result, after all orders in gluon exchanges, is simply

[TABLE]

Lastly, to compute the expectation value over the product of four-dipole operators, we need to sum over each of the elements of the 24-dimensional column vector with the respective $N_{c}$ weights of each of the $n$-tupoles. The four-dipole configuration in this formulation has weight unity. The dipole-dipole-quadrupole configuration comes with $1/N_{c}$, both quadrupole-quadrupole and dipole-sextupole configurations with $1/N_{c}^{2}$, and octupoles with $1/N_{c}^{3}$. The sum over the 24-dimensional basis can then be written as the scalar product of the corresponding row vector and the column vector $\vec{\mathcal{N}}$:

[TABLE]

Here the $\cdots$ denote the different permutations of each of the five configurations in Fig. 4. With this expression for $\mathcal{N}$, the combined gluon exchange and tadpole contributions in Eq. (17) can be written as

[TABLE]

As we noted, analytical expressions for these quantities are too cumbersome to compute. However, with the procedure outlined, they can be computed numerically and utilized to compute the four-particle correlation functions we shall discuss further in the next section.

III.4 Abelian limit

We will consider here the computation of the $m$-dipole expectation value in an Abelianized version of the MV model. In the $U(1)$ theory, the Wilson line again represents the multiple scattering of a charged particle off a classical field Bjorken:1970ah :

[TABLE]

However, here the Wilson line is a scalar valued function, not an $SU(N_{c})$ valued matrix as in the non-Abelian case; this simplifies computations enormously. Expectation values of multiple Wilson lines can be evaluated analogously to the non-Abelian case by first expanding each exponential to second order in the gauge field, followed by replacing all pairwise contractions of the gauge field with the Gaussian expectation value, as in Eq. (13):

[TABLE]

The dipole expectation value is then

[TABLE]

and the correlator of two dipoles can be expressed as

[TABLE]

Similarly, for four dipoles, one obtains

[TABLE]

Higher-point correlators can be found similarly by summing over all combinations of two-point functions; a negative sign is introduced for combinations of two quarks or two anti-quarks.

IV Results

We will now discuss the calculation of $m$-particle production using the expression in Eq. (9). The inputs into this expression include the parameter $B_{p}$ representing the transverse area of the projectile and the function $L(\mathbf{x}_{\perp},\mathbf{y}_{\perp})$ that represents the correlation of gauge field configurations arising from a Gaussian distribution of color sources. In the MV model, the quantity $L(\mathbf{x}_{\perp},\mathbf{y}_{\perp})$ is related to the expectation value of the dipole operator through the expression

[TABLE]

In this work, we will use the functional form

[TABLE]

For the infrared cutoff in the model, we will use the value $\Lambda=0.241~{}\text{GeV}$. This value is obtained from parametrizations of the dipole amplitude used in dipole model computations of structure functions that are fit to the HERA deep inelastic scattering $ep$ data Albacete:2010sy ; Lappi:2013zma . We have checked that the results are qualitatively unchanged within a physically reasonable range of values of $\Lambda$.

Following Ref. Lappi:2013zma , a model independent saturation scale is defined through the relation

[TABLE]

For the remainder of this work, we will specify values of $Q_{s}^{2}$ rather than $(g^{2}\mu)^{2}$. We should point out that the mapping between $(g^{2}\mu)^{2}$ and $Q_{s}^{2}$ contains an explicit dependence on $C_{F}$, as is transparent from Eq. (35). When we compare results at various $N_{c}$, this scaling with $C_{F}$ is taken into account. The exception is the $U(1)$ case, where we take $C_{F}=4/3$ when relating $Q_{s}^{2}$ to $(g^{2}\mu)^{2}$.

We stress that the saturation scale $Q_{s}$ here is that of the target nucleus. There is no analogous saturation scale of the projectile in the model we are considering. This is a consequence of the simplicity of our treatment of the projectile’s constituents, which are comprised of nearly collinear uncorrelated quarks alone. The corresponding average multiplicity per interaction is unity. This can be seen by explicit integration of Eq. (9) and can be contrasted with the expression for the multiplicity found in $k_{\perp}$-factorization Dumitru:2001ux

[TABLE]

where $Q_{s,p}$ and $Q_{s,T}$ are the saturation scales of the projectile and target respectively. The transverse overlap area $S_{\perp}$ can be identified with the projectile area $B_{p}$ in our model. The absence of a projectile saturation momentum is the biggest shortcoming of the above framework. The equivalent scale controlling the momenta of the incoming quarks is $1/B_{p}$, which is held fixed. While the fact that our model and the experimental data seem to both be independent of the multiplicity (at least qualitatively) may be suggestive of a common physical origin, it is far from clear that this will hold for a more realistic model of the projectile.

In what follows, we will primarily plot quantities as a function of the target saturation scale $Q_{s}$. The target saturation scale is a function of both Bjorken $x$ and the impact parameter. As discussed above, the multiplicity is a logarithmic function of $Q_{s}$. Instead, $Q_{s}$ is better thought of as a proxy for the energy of the collision. In the CGC picture of high energy QCD, $Q_{s}$ grows with the center-of-mass energy.

As a final remark, we stress that we only expect to make a qualitative comparison with data. In addition to the shortcomings of the model addressed above, our final state distributions are for quarks and not for hadrons. Any correlations computed in this model will be reduced through a variety of effects, such as fragmentation, and quantum evolution of parton distributions. They therefore provide an upper limit for azimuthal correlations in initial state frameworks.

IV.1 Multiparticle azimuthal cumulants and harmonics

Multiparticle correlations carry a wealth of information on the dynamics of the colliding system. For reviews, see for instance, Refs. Dusling:2015gta ; Schlichting:2016sqo . Azimuthal correlations, in particular, are sensitive measures of collective dynamics in heavy-ion collisions. For systems undergoing collective flow, one can define the $n$th-Fourier harmonic coefficients

[TABLE]

where $\phi$ is the azimuthal angle of a produced particle and $\Psi_{R}$ is the angle of a reaction plane determined by the collective geometry of particles produced in the collision. The determination of a suitable reaction plane may be sensitive to a variety of so-called non-flow contributions such as resonance decays, to give one example. The effect of these non-flow correlations can be minimized using cumulant method Borghini:2001vi , which is now widely used in experimental studies of multiparticle correlations. Two- and four-particle cumulants are defined as

[TABLE]

where the average $\langle\cdots\rangle$ in $c_{n}\{m\}$ is the average of all possible combinations of $m$ particles in an event, followed by an averaging over all events Borghini:2001vi ; Dumitru:2014yza ; Khachatryan:2015waa .

The cumulants above can be expressed in terms of the $m$-particle inclusive distribution by first defining a quantity $\kappa_{n}\{m\}$ as

[TABLE]

corresponding to the $n^{th}$ harmonic of the $m$-particle distribution. The averages can be written as a ratio of the $n^{th}$ azimuthal angle moment of the $m$-particle inclusive distribution normalized by the zeroth moment (the $m$-particle inclusive multiplicity):

[TABLE]

Fourier coefficients are then defined from the above cumulants as

[TABLE]

The motivation for the above definitions becomes transparent under the assumption that the $m$-particle distribution factorizes into a product of single-inclusive distributions correlated with each other only through the event plane angle. This is indeed what one would expect if the system undergoes hydrodynamic flow. In such a framework, the two-particle Fourier harmonic $v_{n}\{2\}$ would, as above, be the square root of the corresponding two-particle azimuthal angle cumulant, the four-particle Fourier harmonic $v_{n}\{4\}$, the fourth root of the four particle azimuthal angle cumulant, and so on. The negative sign in the latter case is appropriate because one anticipates intrinsic four-particle angular correlations to be subdominant relative to the square of the two-particle cumulant in Eq. (41). While the observation of $m$-particle Fourier harmonics is suggestive of some form of collective behavior, we will argue instead that it can result from any physical mechanism where higher cumulants are suppressed relative to the mean and variance of the distribution.

In Fig. (7), we show $v_{n}\{2\}$ as a function of $Q_{s}^{2}$ for Fourier harmonics $n=2,3,4,5$. For each harmonic, we studied the sensitivity of the result for multiple values of the maximum integrated transverse momentum ($p_{\perp}^{\rm max}$). We observe that even for $p_{\perp}^{\rm max}=2$ GeV, the result is insensitive to the high momentum cutoff. Due to the small $x$ evolution of partons in the target, its saturation momentum $Q_{s}$ will increase with the increasing center-of-mass energy of the collision. Since this is the only energy dependent variable in our framework, our result indicates that the $v_{n}\{2\}$ are only weakly dependent on the energy. We also observe that the $v_{n}\{2\}$ have a clear hierarchy with $n$, similar to what is seen in experiment Aad:2014lta .

The four-particle cumulant $c_{2}\{4\}$ is shown as a function of $Q_{s}^{2}$ in Fig. (7). We clearly see that by $Q_{s}^{2}=0.3~{}\text{GeV}^{2}$ that there is a change in the sign of the signal, from positive to negative values, resulting in a real $v_{2}\{4\}=(-c_{2}\{4\})^{1/4}$. The magnitude of the signal only weakly depends on the maximum integrated momentum $p_{\perp}^{\rm max}$. The relatively weak variation in the signal above $Q_{s}^{2}\sim 1$ GeV2 is in qualitative agreement with the experimental findings on the center-of-mass energy dependence of the four-particle cumulant Chatrchyan:2013nka ; Abelev:2014mda ; Khachatryan:2015waa ; Belmont:2017lse ; Aaboud:2017acw .

As we noted previously, a positive definite value for $v_{2}\{4\}$ is natural in hydrodynamic models. However, we know of only two 3+1-dimensional numerical hydrodynamic computations of $v_{2}\{4\}$ in $pA$ collisions Bozek:2012gr ; Kozlov:2014fqa . In the case of Ref. Kozlov:2014fqa , while some of the systematics of $pA$ collisions is reproduced, the model is unable to reproduce similar event-by-event systematics of flow in $AA$ collisions.

There also exist qualitative arguments for positive definite $v_{2}\{4\}$ Dumitru:2014yza in an initial state “color domain” model Kovner:2010xk ; Dumitru:2014dra ; Dumitru:2014vka ; Dumitru:2015cfa ; Skokov:2014tka . Though this model provides an intuitive picture of how multiparticle azimuthal cumulants may be generated in an initial state framework, it is unclear how the oriented background color-electric fields are created from first principles Lappi:2015vta .

In Fig. (8), we plot the ratio of the four- to two-particle integrated $v_{2}\{m\}$. For both values of $p_{\perp}^{max}$, this ratio rises with $Q_{s}^{2}$ and saturates above $Q_{s}^{2}\sim 1$ GeV2. The values obtained $v_{2}\{4\}/v_{2}\{2\}\sim 0.5$–$0.6$ are remarkably close to those measured by the CMS Collaboration that shows this quantity increasing with centrality from $\sim 0.675$ to $0.775$ Chatrchyan:2013nka ; Giacalone:2017uqx . A detailed study of this ratio, and other such ratios, and their centrality dependence, in various models of initial state spatial eccentricities, can be found in Ref. Giacalone:2017uqx .

We now study the multiparticle cumulants differentially in transverse momentum. If we keep the transverse momentum of one particle in Eq. (42) fixed and integrate over the momenta of the remaining $m-1$ particles, the two- and four-particle differential Fourier harmonics are defined as Chatrchyan:2013nka

[TABLE]

where $d_{n}\{m\}$ are the differential analogs of $c_{n}\{m\}$. Figure 9 shows $v_{2}\{2\}$ (left) and $v_{2}\{4\}$ (right) as a function of $p_{\perp}$ for two representative values of the saturation scale $Q_{s}^{2}=1$ and $2$ GeV2. Note that while the results are differential in $p_{\perp}$ for one of the particles the remaining particles are integrated up to a $p_{\perp}^{\rm max}=3$ GeV. We note that for $Q_{s}^{2}=2~{}\text{GeV}^{2}$ our results exhibit behavior similar to that seen in the experimental $p$-Pb data data Chatrchyan:2013nka ; Aad:2013fja . We know of one other theory computation of $v_{2}\{4\}(p_{\perp})$ in small systems Iancu:2017fzn .

Event-by-event fluctuations of $v_{n}$ with $v_{n^{\prime}}$ for $n\neq n^{\prime}$ can be captured by symmetric cumulants Bilandzic:2013kga defined as

[TABLE]

These mixed harmonic angular expectation values are defined analogously to those in Eq. (42) and Eq. (43). They were originally introduced in hydrodynamic frameworks as a measure of the nonlinear response of the initial geometry of the system Bilandzic:2013kga . These have been studied mainly in the context of heavy-ion collisions ALICE:2016kpq . In these heavy-ion systems, the symmetric cumulants describe how correlations between the initial moments of the eccentricities, which are driven by the overlap geometry and thus the centrality of the collision, are carried to the final state. An example of symmetric cumulants are the correlations between the $v_{2}$ and $v_{3}$ azimuthal harmonics, denoted by $\text{SC}(2,3)$. From geometrical considerations, there should be an anti-correlation between the initial ellipticity and the triangularity. When converted to correlations of final state momentum anisotropies, this results in a negative $\text{SC}(2,3)$. Studies of the symmetric cumulants $\text{SC}(2,3)$ and $\text{SC}(2,4)$ for heavy-ion collisions have been carried out in hydrodynamic simulations Giacalone:2016afq ; Zhu:2016puf and in a hadronic transport model Zhou:2015eya .

However, in small systems, this picture should not hold. The initial eccentricities are not believed to be strongly centrality driven, but instead are most likely in response to sub-nucleonic fluctuations Welsh:2016siu . We will touch on the topic of sub-nucleon fluctuations in the next section. In our model, we include sub-nucleon scale fluctuations in sampling the positions of the quarks in the projectile given through our Gaussian Wigner function, Eq. (4). We reported in Dusling:2017dqg that our model produces the correct sign for $\text{SC}(2,3)$ and $\text{SC}(2,4)$ and is in qualitative agreement with the data CMS:2017saf . In Fig. (10), we show predictions for the behavior of $\text{SC}(2,5)$, $\text{SC}(3,4)$, and $\text{SC}(3,5)$. Results for these symmetric cumulants, for heavy-ion collisions alone, were shown at the Quark Matter 2017 conference YZhou:QMtalk . Our results for these cumulants agree qualitatively the results presented. We are unaware of any other theory predictions for these cumulants in light-heavy ion collisions. While not their designated purpose, symmetric cumulants in small systems may be an effective way to access information on initial state sub-nucleon fluctuations.

V Detailed systematics of the results

V.1 Role of the projectile

It is expected that sub-nucleon scale fluctuations play an important role in small systems; hydrodynamic computations including such fluctuations have been performed for $pA$ collisions Mantysaari:2017cni . Thus it is also interesting to ask what similarities our model bears to a constituent quark model based picture. To mock up this effect, we introduce a hard distance cutoff (either minimum or maximum) between all quarks in the amplitude and similarly between all anti-quarks in the complex-conjugate amplitude. This is in addition to the Gaussian sampling of the quark positions from the Wigner function introduced in Eq. (4).

The effect of such a cutoff on $v_{2}\{2\}$ is shown in Fig. (11). Starting with the standard Wigner function in Eq. (4) as a reference, we see that by introducing a minimum distance criteria (separating the quarks by at least a distance of $d_{min}=B_{p}/8$ or $B_{p}/4$) the correlations decrease. This is to be expected because forcing the transverse positions of the quarks to be farther away from each other ensures that they are less likely to interact with the same color domain in the target. We would then expect, on the same grounds, that if we confined the quarks to be at most a distance $d_{max}=B_{p}/4$ or $B_{p}/2$ apart, we would see an increase in the strength of the correlation. This expectation is confirmed by the results shown in Fig. (11).

V.2 $Q_{s}^{2}B_{p}$ dependence

In our model, the parameter $B_{p}$ controls the mean transverse area of the projectile, and therefore the transverse overlap area of the scattering off the target. The scale $1/Q_{s}^{2}$ sets the scale for the size of color domains in the target. Therefore the dimensionless product $Q_{s}^{2}B_{p}$ can be interpreted as the number of domains in the target that overlap the projectile. In fact, this dimensionless parameter controls the strength of all correlations. Fig. (12) shows $c_{2}\{2\}$ for three values of $B_{p}$ for $p_{\perp}^{\rm max}=3~{}\text{GeV}$; note that $B_{p}=4~{}\text{GeV}^{-2}$ is used elsewhere in this work. The inset in Fig. (12) shows the same quantity plotted as a function of the dimensionless scale $Q_{s}^{2}B_{p}$ demonstrating that all results fall onto a universal curve, as they must.

One might expect that the strength of the correlation should fall with the number of domains $1/(Q_{s}^{2}B_{p})$, which is a feature of independent cluster models. However, a falloff that goes like $1/(Q_{s}^{2}B_{p})$ is not seen in the results presented above. The reason is that there is another scale, $p_{\perp}^{\rm max}$, which controls the maximal momentum kick from the target to the probe. The inverse of $p_{\perp}^{\rm max}$ is therefore the smallest distance in the target resolved by the probe. One can therefore construct two dimensionless combinations $Q_{s}^{2}B_{p}$ and $Q_{s}^{2}/({p_{\perp}^{\rm max}})^{2}$; the dependence of our results on the number of color domains also depends on what $Q_{s}^{2}/({p_{\perp}^{\rm max}})^{2}$ is.

For $({p_{\perp}^{\rm max}})^{2}\gtrsim Q_{s}^{2}$, the probe resolves a transverse area in the target that is on the order or smaller than the size of a color domain. Because the probe can resolve the structure within individual domains, we expect to see a falloff in correlations to go approximately as $1/(Q_{s}^{2}B_{p})$. Fig. (13) shows $c_{2}\{2\}$ for $p_{\perp}^{\rm max}=10,20,40$ GeV, all of which satisfy the scaling form $(Q_{s}^{2}B_{p})^{-0.95}$ at large $Q_{s}^{2}$. On the other hand, for $({p_{\perp}^{\rm max}})^{2}\leq Q_{s}^{2}$, the probe only resolves transverse sizes larger than the typical domain size. For these smaller values of $p_{\perp}^{\rm max}$, increasing $Q_{s}^{2}B_{p}$, and therefore the number of color domains, does not change the signal since the probe cannot resolve the change in the number of color domains. Fig. (13) shows that for $p_{\perp}^{\rm max}=3,5$ GeV we see a rather modest falloff with the number of domains $(Q_{s}^{2}B_{p})^{-0.18}$.

Our results in Fig. (13) suggest more generally that for small values of $p_{\perp}^{\rm max}$ relative to the saturation scale $Q_{s}$ azimuthal cumulants in initial state models are weakly dependent on the number of clusters. This independence of Fourier harmonics on the number of clusters has been interpreted previously as occurring due to the collective response of the system Basar:2013hea . While coherent multiple scattering may be collective, it is not a final state effect in $pA$ collisions; the interaction with the target is instantaneous and the scattered quarks do not subsequently rescatter.

V.3 $N_{c}$ dependence

In Ref. Dusling:2017dqg , we showed that an Abelianized version of our model demonstrates the systematics often attributed to collective behavior, $v_{2}\{4\}\approx v_{2}\{6\}\approx v_{2}\{8\}$. Given this finding, it is natural (and of intrinsic interest) to determine the $N_{c}$ dependence of the two- and four-particle azimuthal correlations. The dependence of $v_{2}\{2\}(p_{\perp})$ for $N_{c}=2,3$ was discussed previously in Ref. Lappi:2015vta . In the left panel of Fig. (14), we plot the dependence of integrated $v_{2}\{2\}$ (up to a $p_{\perp}^{\rm max}=2$ GeV) as a function of $Q_{s}^{2}$ for the Abelian ($N_{c}=1$) case and for $N_{c}=2-5$. For large $Q_{s}^{2}$, we observe a convergence of the results for $N_{c}\geq 3$. In the right panel of Fig. (14), we plot the $N_{c}$ dependence of $v_{2}\{4\}$ as a function of $Q_{s}^{2}$. When $v_{2}\{4\}$ is real and large, we expect the second term in Eq. (41) for the four-particle cumulant to dominate. This should then give $v_{2}\{4\}\sim 1/{N_{c}}$. We see from Fig. (14) that $N_{c}\,v_{2}\{4\}$ begins to converge for $N_{c}\geq 3$; however, due to limited statistics, the error bars are large.

We previously reported in Ref. Dusling:2017dqg on results for the symmetric cumulants for $\text{SC}(2,3)$ and $\text{SC}(2,4)$ which were in qualitative agreement with experimental results CMS:2017saf . In Fig. (15), we show the $N_{c}$ dependence of the symmetric cumulants. We find that these symmetric cumulants are extremely sensitive to $N_{c}$. (We have chosen here $p_{\perp}^{\rm max}=2$ GeV.) For the Abelian case, the result is an order of magnitude larger than the finite $N_{c}$ results. Further, $\text{SC}(2,3)$ is positive in the Abelian case, which is not observed in any experimental observations. Interestingly, within the limited number of observables that we have studied, this appears to be the only place where the Abelian version of our model qualitatively differs from the non-Abelian results.

One may infer that there is something specific to the non-Abelian nature of coherent scattering that drives $\text{SC}(2,3)$ to become negative. Going to $N_{c}=2$, we find that $\text{SC}(2,3)$ is identically zero. This is not surprising, as for $N_{c}=2$, all odd harmonics are identically zero. This is analogous to the absence of odd harmonics for gluons scattering off a target Schenke:2015aqa . The underlying reason is that $SU(2)$ representations are real, regardless of whether they are in the fundamental or the in adjoint. For $SU(3)$, only the adjoint representation is real. Thus one expects qualitatively different results for even-odd cumulants for $N_{c}\geq 3$ relative to the Abelian model and for $N_{c}=2$.

V.4 Comparison to Glasma graphs

To elucidate the mechanism responsible for the observed negativity of $c_{2}\{4\}$ for our model, we compare this result to that from the Glasma graph approximation. In this Glasma graph approximation, which is applicable for $p_{\perp}>Q_{s}$, the Wilson lines are expanded out to lowest nontrivial order in the gauge fields. The Glasma graph approximation is reviewed in the Appendix, where we also compute the four-dipole correlation function in this approximation. This approximation was used previously to successfully describe two-particle correlations Dumitru:2008wn ; Dumitru:2010iy ; Dusling:2012iga ; Dusling:2012cg ; Dusling:2012wy ; Dusling:2013qoz ; Dusling:2015rja , especially near side “ridge” correlations.

Fig. (16) shows a comparison of $c_{2}\{4\}$ in the Glasma graph approximation to our result, which includes all order contributions from the Wilson lines. It is clear that coherent multiple scattering in the MV model computation for $c_{2}\{4\}$ drives it to become negative. In contrast, the Glasma graph approximation to this full result is always positive.

It is interesting to explore further the physics underlying this striking result. Intrinsic $n$-particle correlations in the Glasma graph approximation are large. Indeed, the strength of the correlated piece in this distribution relative to the disconnected product of $n$ single-particle distributions is the same for all $n$; this is close to that of an $n$-particle Bose distribution Gelis:2009wh . Our results suggest that coherent multiple scattering depletes these higher-point intrinsic correlations. In Eq. (41) for instance, this would lead to the second term (the square of $c_{2}\{2\}$) to dominating over the first term, which comes from intrinsic four-particle correlations.

V.5 Rapidity dependence

Before we conclude our discussion of features of the model, it is appropriate to discuss how rapidity correlations manifest themselves in this framework. Since long range rapidity correlations are an essential feature of the experimentally observed ridge-like correlations, it is important to determine whether such correlations are present in this framework. This is particularly so since the hybrid formalism employed in analytical studies of such multiparticle correlations Dumitru:2002qt ; Kovner:2012jm ; Kovchegov:2012nd is valid in the forward rapidity region. More precisely, $x$ in the projectile should be relatively large, with typical values for large-$x$ taken to be $x\geq 0.01$. However this does not imply that the resulting correlations are short-range in rapidity. We will show this explicitly by reintroducing rapidity dependence into the single-particle and multi-particle distributions.

We consider an eikonal quark in the projectile traveling in the $z^{+}$ direction with initial momentum $k^{\mu}=(k^{+},0,{\bf k_{\perp}}={\bf 0})$ and final momentum $p^{\mu}=(p^{+},0,{\bf p_{\perp}})$ after its scattering off the target. In the hybrid framework, the differential multiplicity of the scattered quark with the final state of momentum ${\bf p}\equiv(p^{+},{\bf p_{\perp}}$) is given as

[TABLE]

The above result can be worked out following the formalism of Ref. Dumitru:2002qt . (Note though that in Ref. Dumitru:2002qt , the quark is traveling in the $z^{-}$ direction.) The expression for $dN/d^{2}{\bf p_{\perp}}$ is given by Eq. (6), where averaging over the target we previously defined is implicitly assumed. Therefore, for single-quark scattering the distribution is a delta function in $\delta(p^{+}-k^{+})$, which is simply a consequence of the eikonal approximation.

The single inclusive distribution of quarks produced in $pA$ collisions is obtained by convoluting the above expression with the quark parton distribution, which represents the probability of finding a quark in the proton wavefunction:

[TABLE]

To obtain the single inclusive distribution of hadrons, this expression has to further convoluted with a fragmentation function. This will quantitatively modify the rapidity dependence but will not modify it qualitatively. We will therefore not further consider this point.

The longitudinal momentum carried by the initial quark is $k^{+}=\frac{\sqrt{s}}{\sqrt{2}}x_{q}$, where $x_{q}$ is the quark’s momentum fraction. Likewise, the longitudinal momentum of the final-state quark can be written in terms of its rapidity as $p^{+}=\frac{p_{\perp}}{\sqrt{2}}e^{y}$. Substituting Eq. (48) into Eq. (49), we then obtain

[TABLE]

where $x_{q}^{\prime}$ has been set by the $\delta$-function to be

[TABLE]

In Fig. (17(a)), we show the single-particle inclusive quark distribution using the NNPDF NLO singlet PDF at $Q^{2}=9~{}$GeV2 Ball:2014uwa . For $p_{\perp}=3$ GeV and $\sqrt{s}$ = 5.02 TeV, a value of $x_{q}=0.01$ corresponds to a rapidity of $y\simeq 2.8$. Thus, while strictly speaking, our approach is valid for $y\geq 2.8$, we can see that the particle production does extend over a wide range in rapidity and is not peaked in the forward direction.

This can be extended to multiparticle production in the same fashion. For instance, for two particles, we would obtain

[TABLE]

This two-particle distribution is shown in Fig. (17(b)). The rapidity of the first quark is fixed at the edge of where the hybrid approach works. We see that the correlation persists out to a $\Delta y$ of 3 to 4 units as the rapidity of the second quark is varied.

We conclude by pointing out that when computing $v_{n}\{m\}$ we are taking ratios of the momentum integrated $m$-particle distributions, where the numerator is also convoluted with a cosine function. The factorized form of the rapidity dependence in Eq. (52) is then suggestive that the resulting $v_{n}\{m\}$ will be weakly dependent, or perhaps even constant, as a function of rapidity.

VI Conclusions

In this work, we elaborated significantly on a parton model framework for multiparticle correlations that was first presented in Ref. Dusling:2017dqg . In this model, collinear partons that are localized in a transverse area $B_{p}$ in the projectile, scatter off color domains of size $1/Q_{s}^{2}$ in the target. The building blocks in this framework are dipole correlators. For a quark projectile, these correspond to a color trace over a path ordered lightlike Wilson line of target fields at a given transverse spatial position in the amplitude multiplied by its adjoint in the complex-conjugate amplitude, normalized by the number of colors, $N_{c}$. The azimuthal cumulants of $n$-particles are proportional to the expectation value over the product of $n$-dipole correlators. We discussed at length the procedure to compute the expectation value of four-dipole correlators. Gluon exchanges among the dipoles generate distinct quadrupole, sextupole, and octupole topologies and the permutations of their spatial positions, generating a 24 by 24 matrix, that can be exponentiated to determine the expectation value of four-dipole correlators. Our procedure can be extended straightforwardly to compute expectation values of products of an arbitrary number of dipole operators.

We presented results for $v_{n}\{2\}$ as a function of $Q_{s}^{2}$, demonstrating a clear hierarchy in the $n=2,3,4,5$ harmonics. These results are only weakly dependent on $p_{\perp}^{\rm max}$ upper limit in the integrals. All the harmonics show only a weak dependence on $Q_{s}^{2}$. Since $Q_{s}^{2}$ in the CGC framework increases with increasing energy, and centrality, our results are only weakly dependent on these. Note further that our results, by construction, are independent of the multiplicity. These results for $v_{n}\{2\}$ exhibit the qualitative features of the data seen in light-heavy ion collisions at RHIC and the LHC.

We next presented results for $c_{2}\{4\}$ as a function of $Q_{s}^{2}$ for two different values of $p_{\perp}^{\rm max}$. For both these values, $c_{2}\{4\}$ changes sign around $Q_{s}^{2}=0.3$ GeV2 and becomes increasingly negative before appearing to saturate. For large $Q_{s}^{2}$, the results are sensitive to $p_{\perp}^{\rm max}$. The negative value of $c_{2}\{4\}$ corresponds to a real $v_{2}\{4\}$. We computed the ratio of $v_{2}\{4\}/v_{2}\{2\}$. In hydrodynamic models, such ratios are sensitive to the initial geometry in the system, motivating experimental extractions of the same. The values we obtained are about 10% lower than the data, which at present have significant error bars.

The dependence of $v_{2}\{2\}$ and $v_{2}\{4\}$ as a function of $p_{\perp}$ increases with $p_{\perp}$ before saturating and turning over. For both quantities, this saturation occurs later with increasing $Q_{s}^{2}$. In particular, for $v_{2}\{4\}$, we observed for $Q_{s}^{2}=2$ GeV2 that it is quite flat for the $p_{\perp}$ range between 1 and 2 GeV. These features of our results are also qualitatively similar to data on small-particle systems. To the best of our knowledge, no computations exist in other models for $v_{2}\{4\}(p_{\perp})$.

Symmetric cumulants $\text{SC}(n,n^{\prime})$, which measure the correlation of $n$th Fourier harmonics with $n^{\prime}$ Fourier harmonics, were constructed to understand the nonlinear hydrodynamic response of the system to correlations in the initial spatial geometry. We studied these in our initial state framework as a function of $Q_{s}^{2}$. We showed in Ref. Dusling:2017dqg that $\text{SC}(2,3)$ and $\text{SC}(2,4)$ computed are in qualitative agreement with the data presented for heavy-ion collisions and in light-heavy ion collisions. Here, we make predictions for the $\text{SC}(3,4)$, $\text{SC}(2,5)$, and $\text{SC}(3,5)$ cumulants.

We examined closely the dependence of our results on $Q_{s}^{2}B_{p}$ and $Q_{s}^{2}/({p_{\perp}^{\rm max}})^{2}$, the two dimensionless parameters in our model. The former corresponds to the number of color domains in the target that are encountered by the projectile. The latter corresponds to the resolution of the partons in the projectile to the structure of the color domains. Interestingly, we find that for larger values of $Q_{s}^{2}/({p_{\perp}^{\rm max}})^{2}$, the two-particle cumulants are only weakly dependent on the number of color domains. In contrast, for smaller values of $Q_{s}^{2}/({p_{\perp}^{\rm max}})^{2}$, we find that the cumulant behaves approximately as $1/(Q_{s}^{2}B_{p})$, as would be anticipated in an independent cluster model. Our results suggest therefore that the $p_{\perp}^{\rm max}$ considered is important in any interpretation of the data that may be construed as satisfying or violating an independent cluster model.

We studied next the dependence of our results on $N_{c}$. In Dusling:2017dqg , we showed that the Abelianized treatment of our model reproduced the pattern of $v_{2}\{2\}>v_{2}\{4\}\approx v_{2}\{6\}\approx v_{2}\{8\}$ seen in the data on $pA$ collisions at the LHC. We studied further the $N_{c}$ dependence of $v_{2}\{2\}$ and $v_{2}\{4\}$ and demonstrated that the both behave as $1/N_{c}$ for $N_{c}\geq 3$. There is therefore no ordering in $N_{c}$ among $m$-particle cumulants. Practically, it means that $N_{c}$ suppressed topologies in products of lightlike Wilson lines must be included in such computations.

Unlike the azimuthal cumulants $c_{2}\{4\}$, the symmetric cumulants $\text{SC}(n,n^{\prime})$ (where $n$ and $n^{\prime}$ denote distinct Fourier harmonics) have a qualitatively different behavior in the Abelian formulation of the model relative to those for the $N_{c}=3$ case we studied previously in Ref. Dusling:2017dqg . This qualitative difference is not unique to $N_{c}=1$. It is also seen for $N_{c}=2$. In this latter case, the odd harmonics are strictly zero; hence, the corresponding symmetric cumulants are also zero. As we discussed, the underlying reason is that for $SU(2)$ both fundamental and adjoint representations are real. This is responsible for two-particle correlations being symmetric about relative azimuthal angles of zero and $\pi$.

To obtain deeper insight into our results, we examined our results within the Glasma graph approximation to our framework. This approximation corresponds to expanding out the Wilson lines in the dipole correlators to lowest nontrivial order. Physically, it corresponds to each quark line interacting at most with two gluons, either in the amplitude, or the complex-conjugate amplitude, or across the cut separating the two. It is justified when $Q_{s}^{2}\ll p_{\perp}^{2}$. Remarkably, we find that our results for $c_{2}\{4\}$ are qualitatively different when we include coherent multiple scattering ($Q_{s}^{2}/p_{\perp}^{2}$ to all orders) as opposed to the Glasma graph approximation. In the former case, one obtains a real and positive $v_{2}\{4\}$; in the latter case, one does not. Therefore, the absence of $v_{n}\{m\}$ multiparticle correlations in previous Glasma graph treatments is an artifact of the approximation and not a genuine feature of initial state frameworks.

We noted that multiparticle correlations are quite strong in the Glasma graph approximation, and are close to those of an $m$-particle Bose distribution. Indeed, this may be the reason why one does not see signatures of “collectivity”, as defined by the $v_{n}\{m\}$ Fourier harmonics. These assume that the mean and variance dominate the cumulant distribution. Our results suggest that coherent multiple scattering dilutes the contributions from the higher cumulants relative to the mean and the variance, thereby generating the aforementioned signatures of collectivity. In our model, however, the coherent scattering is virtually instantaneous. It takes place on a time scale corresponding to the time it takes partons to cross a Lorentz contracted nucleus. Further, the partons that scatter off the common color field do not rescatter.

The origin of this putative signature of collectivity therefore has little to do with hydrodynamics per se. However, our results do not exclude the possible presence of final state interactions, or even hydrodynamics, in the data on light-heavy ion systems. They do provide a clear and simple counterexample to interpretations of these signatures as being of unique origin. Such signatures of collectivity must also be consistent with other signatures of collectivity. In the larger heavy-ion collision systems, jet quenching is seen very clearly, and is an independently robust measure of final state interactions.

It is interesting to consider how this model can be developed further. Gluon degrees of freedom, which are not apparent at lower energies, become important when hadrons are boosted to higher energies. As is well known, a bremsstrahlung cascade develops, which then has a shock wave interaction with the nucleus. The partons in the cascade subsequently fragment to hadrons. This picture is implicit in the CGC+PYTHIA model of so-called dense-dense collisions of small systems Schenke:2016lrs that combines Yang-Mills dynamics of gluons Schenke:2015aqa with Lund fragmentation Andersson:1983ia . Because event-by-event simulations are essential to compute multiparticle cumulants, such computations are computationally intensive. Our work suggests the need for further developments in this direction.

VII Acknowledgements

We would like to thank Jiangyong Jia, Tuomas Lappi, Aleksas Mazeliauskas, Larry McLerran, Jean-Yves Ollitrault, Jean-François Paquet, Björn Schenke, Sören Schlichting, Juergen Schukraft, Chun Shen, Vladimir Skokov, and Prithwish Tribedy for useful discussions. This material is based on work supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under Contracts No. DE-SC0012704 (M.M. and R.V.) and No. DE-FG02-88ER40388 (M.M.). M.M. would also like to thank the BEST Collaboration for support. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231 and the LIRED computing system at the Institute for Advanced Computational Science at Stony Brook University.

Appendix A Appendix: Glasma graph approximation

In this section, we show how our present calculation fundamentally differs from the so-called “Glasma graph” result Dumitru:2008wn ; Dumitru:2010iy ; Dusling:2012iga ; Dusling:2012cg ; Dusling:2012wy ; Dusling:2013qoz ; Dusling:2015rja ; Ozonder:2014sra ; Ozonder:2016xqn . The Glasma graph approximation is constructed by considering all possible two-gluon exchanges between quarks comprising the projectile and the target nucleus under the assumption of Gaussian statistics. The expectation value of $n$-Wilson lines in the Glasma graph approximation can be evaluated by expanding the path ordered exponential of the Wilson line to order $n$ in the coupling constant. The resulting expectation values of gauge fields are then evaluated by re-expressing higher order expectation values as a product of two-point functions. This procedure was followed in Ref. Lappi:2015vta in order to evaluate the dipole-dipole correlator in the Glasma graph approximation. Here, we extend the derivation and results of Ref. Lappi:2015vta to higher order correlators but take a more diagrammatic approach.

Diagrammatically, the Glasma graph approximation amounts to replacing each two-gluon exchange with the expectation value of a single dipole operator. For example, for single-quark scattering we find the relation of Fig. (18).

While the tadpole terms are explicitly shown in the single scattering case in Fig. (18), for higher-point functions diagrams containing two gluons on the same quark are power suppressed by either $k_{\perp}/Q_{s}\ll 1$ or $1/\left(BQ_{s}^{2}\right)\ll 1$. (A line containing two quarks on opposite sides of the cut–corresponding to a gluon connecting a quark with its conjugate amplitude at a different coordinate–is not suppressed and therefore included in what follows.) While these additional terms are formally higher order and will be ignored in the discussion to follow, we should point out that they were found to be important in obtaining a quantitative agreement with experiment Dusling:2012cg .

One can generate all possible Glasma graphs by starting with the completely disconnected diagram whereby each quark multiple scatters independently off the target. We show an example of this completely disconnected contribution for two-quark scattering as the leftmost diagram of Fig. (19).

Further diagrams are generated by finding unique exchanges of coordinates (exchanging two gluon end points) along with an accompanying $\left(N_{c}^{2}-1\right)$ suppression for each exchange. There are two unique topologies for two quark scattering as shown in Fig. (19). The resulting expression for the Glasma graph approximation to the expectation value of two-dipole operators is therefore

[TABLE]

Higher-point correlators can be found in a similar fashion. For $n$ quarks there are $(n-1)!!\equiv n(n-2)\cdots 3\cdot 1$ diagrams resulting from all possible unique contractions between quark lines. For a correlation among six Wilson lines there are $5\cdot 3\cdot 1=15$ pair-wise contractions while for a correlator among eight Wilson lines there are $7\cdot 5\cdot 3\cdot 1=105$ pair wise contractions.

As an aside, we remind the reader that the combinatorics discussed above are relevant for a dilute-dense framework. A Glasma graph approximation has also been employed in the dense-dense limit which was shown to have $(n-1)!!^{2}$ diagrams for $n$-gluon production. In the dense-dense limit there would be $15^{2}=225$ and $105^{2}=11025$ diagrams for the six- and eight-point functions respectively as shown previously in Refs. Dusling:2009ar and Gelis:2009wh .

We now come to the expectation value of four dipoles in the Glasma graph calculation following the procedure outlined above. The starting point is the completely disconnected contribution in which each of the four quarks scatters independently as shown in the left diagram in Fig. (20).

Starting from the disconnected diagram, there are 12 unique coordinate exchanges that can be made. One example is shown in the right diagram of Fig. (20) where the coordinates $\bar{\mathbf{x}}$ and $\bar{\mathbf{w}}$ have been swapped. These 12 diagrams resulting from a single coordinate exchange result in a contribution of $\left(N_{c}^{2}-1\right)^{-1}\mathcal{T}_{1}$ where

[TABLE]

We next consider diagrams resulting from unique and non-trivial two-coordinate exchanges. There are two classes of diagrams which enter at this order. The first is shown in the left diagram of Fig. (21). The diagram factorizes into two sub-graphs, one being a completely connected three-quark scattering and the other being a single independent quark scattering. There are 32 such diagrams which will contribute with a factor of $\left(N_{c}^{2}-1\right)^{-2}\mathcal{T}_{2a}$ where

[TABLE]

The second class of diagrams is shown in the right of Fig. (21). It again factorizes into two sub-graphs, but in this case, each sub-graph is a completely connected two-quark scattering event (i.e. one of the two connected graphs shown in Fig. (19)). There are 12 such diagrams which will contribute as $\left(N_{c}^{2}-1\right)^{-2}\mathcal{T}_{2b}$ where

[TABLE]

The last set of diagrams corresponds to those obtained by three unique and non-trivial coordinate exchanges. An example of such a diagram is shown in Fig. (22). There are 48 such diagrams which contribute as $\left(N_{c}^{2}-1\right)^{-3}\mathcal{T}_{3}$ where

[TABLE]

Putting all of the diagrams together, we arrive at the final expression for the expectation value of four dipoles in the Glasma graph approximation,

[TABLE]

We now come to the evaluation of the four-particle cumulant using the above Glasma graph approximation in the four-particle inclusive distribution, where the four-particle cumulant is defined in Eq. (41). The Glasma graph approximation implicitly assumes that $\left(BQ_{s}^{2}\right)\gg 1$, and one must take this power counting into consideration when taking the ratios in the above expression for the cumulants. In order to see this more clearly let us start with the expression for the double inclusive distribution,

[TABLE]

We can evaluate the single inclusive multiplicity by integrating over $d^{2}p_{1}$ and $d^{2}p_{2}$ resulting in

[TABLE]

In the full nonlinear theory (without approximation), we know that $\left<D\left(R_{1},R_{1}\right)D\left(R_{2},R_{2}\right)\right>=1$ and therefore $N=1$. However, in the Glasma graph approximation, the dipole-dipole correlator is replaced by Eq. (53) resulting in the following expression for the total multiplicity,

[TABLE]

violating unitarity by a term of order $1/\left(BQ_{s}^{2}\right)\ll 1$. When evaluating the cumulant in Eq. (41) one should formally only keep the leading order terms in $1/\left(BQ_{s}^{2}\right)$. In practice, this means retaining only the leading disconnected contribution in the two terms $\kappa_{0}\{4\}$ and $\kappa_{0}\{2\}$ appearing in the denominators; these terms were defined in Eq. (42). For consistency within the Glasma graph approximation, we should therefore take $\kappa_{0}\{4\}=\left(\kappa_{0}\{2\}\right)^{2}=\left(\kappa_{0}\{1\}\right)^{4}$ where $\kappa_{0}\{1\}$ is just the single-inclusive distribution.

Now, coming to the numerator, due to rotational invariance, $\mathcal{T}_{1}$ and $\mathcal{T}_{2a}$ do not contribute to $\kappa_{2}\{4\}$. Out of the 12 diagrams in $\mathcal{T}_{2b}$, four vanish by rotational invariance, and the remaining eight cancel with the term $2\left(\kappa_{2}\{2\}\right/\kappa_{0}\{2\})^{2}$. So, within the Glasma graph approximation, the cumulant can be evaluated using

[TABLE]

where $\kappa_{2}\{4\}$ is evaluated using the 48 diagrams contained within $\mathcal{T}_{3}$. With this, we compare this Glasma graph result with the full non-linear result introduced in Sec. III in Fig. (15) in the main text, with accompanying discussion in Sec. V.4.