General formulae for dipole Wilson line correlators with the Color Glass Condensate

We present general formulae to compute Wilson line correlators with the Color Glass Condensate described by the McLerran-Venugopalan model. We explicitly construct a complete and non-orthogonal set of color-singlet bases and write matrix elements down, so that the exponential of the matrix leads to the Wilson line correlators. We further develop a systematic perturbative expansion of dipole Wilson line correlators in terms of 1/Nc where Nc is the color number. As a phenomenological application we calculate the flow harmonics vn{m} in the dipole model and discuss the Nc scaling.


JHEP11(2017)114
characterized by Q 2 s . The non-linear quantum evolution equation for the color source distribution is known as the JIMWLK equation named after the authors of the pioneering works [12-14, 16, 17]. Although solving the JIMWLK equation demands huge computational resources (see, for example, ref. [18] for the reformulation using the Langevin equation, and refs. [19,20] for numerical simulations), a Gaussian approximated solution has been derived [21,22]. The simplest CGC model in the Gaussian approximation is commonly called the McLerran-Venugopalan (MV) model, which is often used as an initial input for the JIMWLK evolution, or, this model itself could provide us with a good description of the qualitative features of high-energy QCD processes.
One of the most interested and testable quantities calculated in the CGC framework is the particle production [23][24][25][26] and correlation [27][28][29][30][31], including electromagnetic probes [32][33][34][35][36]. In the MV model various correlations among gluons and quarks have been discussed, and these predictions are to be compared to the experimental data of jets and hadron correlations. One of the most well known examples is the CGC picture to understand the ridge structure seen in the di-jet or di-hadron correlations in rapidity space [37,38]. Thus, it would be a natural extension to apply the CGC-based calculations to account for the collective behavior in small size systems. As the collision energy grows up, small size systems such as the p-A (proton-nucleus) and even the p-p collisions may look similar to the A-A collision, yielding systematic patterns of the flow observables (see ref. [39] for an experimental overview). The final state interaction might be still responsible for the flow observables, but at the same time, one must also estimate the initial state effect to quantify which of the initial and the final state interactions is more important. Recently, within the dipole model for the particle production the systematics of the higher order flow observables, v n {m} (i.e. n-th harmonics of the m-particle flow) has been quantified in the glasma graph approximation [40,41] and in the full MV model [42,43].
Generally, in the CGC model calculations, the most time-consuming part in numerics is the numerical computation of the expectation value of the Wilson lines. In the presence of the CGC background fields corresponding to the saturated soft gluons typically in the p-A collision and the forward p-p collision, one needs to take account of multiple scatterings which amount to the Wilson line in the eikonal approximation. It has been repeatedly discussed how to compute the Wilson line correlators efficiently in separate contexts (see ref. [44] for a very recent work toward general formulae for the Wilson line correlators, for instance). Therefore, it would be quite useful to establish a prescription to compute the Wilson line correlators with some generality but in a fairly straightforward way. There are several preceding works along these lines from an early attempt in ref. [45] to a very recent reformulation in ref. [46] on top of an explicit calculation in ref. [44]. The purpose of the present work is to translate the powerful method of ref. [45] to a more specific physics problem in a more handy form. In particular, the setup in ref. [45] was too general and it was not clear how to utilize the formulae for phenomenological applications. In the present work, hence, we focus on the special type of correlation function in terms of the dipole operators in the color fundamental representation, which is the basic building block in the dipole model [47]. For n-point dipole correlations the problem is reduced to the exponentiation of n! × n! matrices. The calculation of the exponential of the matrix JHEP11(2017)114 is numerically feasible but it is sometimes insightful to perform the large-N c expansion, especially to identify the N c scaling of observables. We emphasize that our formulae take a quite convenient form for the large-N c expansion and we will demonstrate the expansion explicitly up to the order where the first nonzero flow cumulants appear. We check the validity of the large-N c formulae by comparing to the exact results known for n = 2 with N c = 3. For n ≥ 4 we need to be very careful of the correct large-N c counting and we will argue the subtlety in the numerical analysis.

Master formulae
The Wilson line correlators that we calculate in this paper are expressed in analogy to quantum mechanics as follows [45]: where the free Hamiltonian, H 0 , is defined as Here, we take the summation over a in group space (which is always implicitly assumed).
In this paper we limit our considerations to the fundamental representation only and T F a i 's represent the elements of su(N c ) algebra in the fundamental representation, while the formula is valid for any representation. Here, L(0, 0) is quadratically divergent, so that only color states that have zero eigenvalues of H 0 can make finite contributions to the Wilson line correlators [and thus, the detailed definition of L(0, 0) is not important here; see ref. [45] for the explicit form of L(0, 0)]. We note that the definition of Q s is slightly different from refs. [42,43] and we will come to this point later when we apply our formulae for the flow observables. For the moment, eq. (2.2) is understood as a definition of Q s in our convention. The color matrix in the interaction part is In the above the common building block, Γ(x ⊥ ), is defined in the CGC formalism by the following integral: where J α (x) is the Bessel function of the first kind, Λ is an infrared (IR) cutoff of order of Λ QCD and we introduced a shorthand notation asΛ := 1 2 Λe γ−1 . This approximated JHEP11(2017)114 expression has undesirable behavior in the IR region especially when |x ⊥ |Λ < 1. We cure this IR problem, according to ref. [43], by introducing another regulator as (2.5) We will adopt this regularized approximation for Γ(x ⊥ ) in our numerical calculations later in section 5. Now we shall sketch our strategy to proceed to concrete calculations of the Wilson line correlators using eq. (2.1). We will first identify all the color-singlet bases with which H 0 vanishes. For the n-th power product of U and U * , there are 2n Wilson lines, and then the number of independent color-singlet bases should be n!, as we will explain in the next section. We can easily construct the color index structures by taking the permutations. After confirming that H 0 surely vanishes with such bases, next, we will consider the matrix elements of V . In general it is not an easy task to find analytical expressions for the eigenvalues of V . Nevertheless, it is known that the large-N c approximation works at good quantitative accuracy, and we will systematically make an expansion in power of 1/N c to find analytical expressions.

Color singlet bases
Because of the dipole-type structures of the Wilson lines it is easy to find all the combinations of the color singlet indices. One trivial singlet is immediately found as {α}; {ᾱ}|s 0 := N −n/2 c δ α 1ᾱ1 · · · δ αnᾱn .
(3.1) Clearly all the permutations are possible, i.e. we introduce |s p as a permutation of |s 0 . For this purpose let us introduce the symmetry group S n with elements π p that denotes a permutation as π p = 1 2 · · · n p 1 p 2 · · · p n . We note that any π p can be expressed as a product of cycles. For example, it is easy to confirm the following relation, where one-cycles, i.e. (3) in the above case, are trivial and not explicitly written. Let us see several useful formulae for later calculation checks. Because of the well-known relation,

JHEP11(2017)114
we can easily prove that which is schematically expressed as For the complex conjugate, because T F a 's are Hermitean, the following should hold: Thanks to the index structure of |s 0 , we readily see; {π p β}; {(i j)β}|s 0 = {π p (i j) −1 β}; {β}|s 0 and trivially (i j) −1 = (i j), so that we can again give a schematic representation as In this case, if p j = i, the contract of δ β iβj δ α iᾱj = δ βp jβ j δ αp jᾱ j in the above and δ αp jᾱ j in |s p makes N c δ βp jβ j . For p j = i the contract of δ β iβj δ α iᾱj and δ αp jᾱ j δ α iᾱ p −1 indicates an index that satisfies Therefore, we establish the following schematic relations: In what follows below, we will utilize the formulae (3.7), (3.10), and (3.12) to perform calculations with a compact notation. It would be an instructive check to see how H 0 |s p = 0 is satisfied for all |s p . For this purpose we first need to expand the second-order Casimir operator in H 0 as

JHEP11(2017)114
Here, again, we note that the summation over a is always implicitly assumed. The first term is nothing but the Casimir operator, so that it is simply given by according to the su(N c ) algebra. Using the formulae (3.7), (3.10), and (3.12), we can simplify the rest (if applied to |s p ) as The first term cancels out with eq. (3.14). Since i and j run from 1 to n, we can equivalently take a summation with respect to j instead of using p j , leading to Also, we see π p · (i j) = (p i p j ) · π p , which allows us to replace p i and p j with i and j in the summations. Finally we arrive at This completes our confirmation of H 0 |s p = 0 for any |s p . For the practical calculation the most important is the evaluation of the matrix elements of V . In our formalism we can infer the matrix elements from which can be easily verified with eqs. (3.7), (3.10), and (3.12). We can further simplify the above expression by introducing a notation for a proper combination of Γ's, i.e., We then define the explicit components of the matrix elements as

JHEP11(2017)114
Here we must be careful of the fact that |s p 's span a complete set of bases but they are not orthogonal. Using these notations and definitions we can summarize the non-zero components as follows: and other matrix elements vanish. These formulae are our central results in the present paper. For the actual application of the formulae, we should compute the exponential of V as seen in eq. (2.1).
To understand how the formulae work, let us consider the simplest example of n = 2. The matrix elements of 2 × 2 matrix V read It is a straightforward calculation to obtain two eigenvalues as λ ± = 1 2 (trV ± ϕ) with Now we can express the exponential of V in a simple form as (3.26) Although the calculation machinery is rather simple, a larger n would cause a huge computational cost. Hence, we will seek for an algorithmic expansion to approximate e −V without complicated matrix algebra.

Dipole Wilson line correlators and the large-N c expansion
For the application for the particle production problem in the relativistic heavy-ion collision [42,43], we are specifically interested in the correlation functions of the dipole operators. The definition of the dipole operator is

JHEP11(2017)114
From this form it is obvious that the n dipole expectation value is given by an matrix element of e −V evaluated with |s 0 , i.e., We may be able to do a direct computation, but we can develop a more sophisticated method assuming that N c is large enough. In view of eqs. (3.21) and (3.22), the offdiagonal components are suppressed by 1/N c , so that we can avoid exponentiating V but make a systematic expansion in terms of V p(ij),p . For notational brevity we shall denote the diagonal and the off-diagonal parts of V as V (0) and V (1) , respectively. Then, the starting point for the systematic perturbative expansion is the interaction picture as in quantum mechanics expressed as where T τ stands for the time-ordered product in terms of τ and the time dependent V (1) (τ ) in the interaction picture is defined as Thus, up to the second order in V (1) for example, the perturbative expansion reads Because V (0) is a diagonal matrix, its matrix elements, V p,p , are the eigenvalues of V . Then, let us introduce an eigenvector |p with an eigenvalue E p for V (0) . That is, where we decomposed the eigenvalue according to the 1/N c order as It is then easy to re-express the first perturbative correction as

JHEP11(2017)114
Here, V qp is defined as V (1) |p = q |q V (1) qp , where we note that this V (1) is an original matrix, not the one in the interaction picture.
In the same way, we can proceed to the second perturbative correction as At this point, we can see a general algorithm to go to arbitrary high orders. The next order, for example, is generated automatically via one more iteration as Now, we are ready to compute s 0 |e −V |s 0 up to the N −2 c order. Noting that V (1) has a matrix element between p and p(ij), we can write down As we already mentioned, |s p 's are not orthogonal for different p's, and a simple calculation leads to the normalization as s 0 |s p = N −n+np c where n p denotes the number of cycles of π p . Because the second and the third terms are already suppressed by N −2 c , we can replace E (ij) with E (0) (ij) in the above expansion. Then, we notice that E (0) (ij) is always accompanied by E (0) 0 , which motivates us to introduce a new notation as (4.13) Now, by expanding E 0 and using the above relations, we can reach the result from the large-N c expansion up to the second order as .

Comparison to the exact answer
In this section let us make a comparison between numerical results from our expansion (4.14) and the exact answer. In particular for the n = 2 case as we discussed around eq. (3.26), the full analytical expression for the dipole Wilson line correlation is known for general N c , which provides us with a useful benchmark to quantify the validity of the large-N c approximation in eq. (4.14). Because our present purpose is to check our formulae (4.14) , concrete values of model parameters are not relevant. The coupling g always appears as a combination of g 4 Q 2 s , so we can take g = 1 without loss of generality and change Q s . We measure all variables in the unit ofΛ here. This means, for example, Q s = 5 in this section is actually Q s = 5Λ, etc. Figure 1 shows the validity check between eqs. (3.26) and (4.14) for the n = 2 dipole correlator with N c = 3. The agreement generally depends on x ⊥i and y ⊥i , but our formulae (4.14) work quite well for almost all x ⊥i and y ⊥i as long as Q s is not too large (when Q s is too large, the outputs are too small, and the errors become relatively larger). Here, in figure 1, we chose Q s = 5 and x ⊥1 = (0, 0), y ⊥1 = (0.6, 0.6), x ⊥2 = (x 2 , 0), y ⊥2 = (0.3, 0.3), which is intentionally chosen to make the difference as visible as possible for the small x 2 region, and so, the overall agreement is better than shown in figure 1.
We next see the Q s dependence of the n = 2 dipole correlator together with the Abelian approximation, as depicted in figure 2. We introduce the Abelian approximation as employed in refs. [42,43] so that the n = 1 expectation value can reproduce the exact result. For example of the n = 2 case, the Abelian approximation reads which would agree with the exact answer in the limit of x ⊥1 = y ⊥1 or x ⊥2 = y ⊥2 (in which the correlator reduces to the n = 1 one), but deviates from the exact answer for general x ⊥i and y ⊥i . In this sense, the expression like eq. (5.1) is to be regarded as an Abelianized JHEP11(2017)114 extrapolation from the n = 1 answer. Figure 2 clearly shows that a small disagreement magnified in figure 1 is actually a negligibly small difference in the whole profile over various Q s . The Abelian approximation captures a qualitative dependence with increasing Q s , while the quantitative values should be considered as only order estimates.

Flow harmonics and higher order contributions
We follow the calculations of the flow observables, v n {m}, in the dilute-dense system according to refs. [42,43], i.e. we focus on only the CGC effects on v n {m} apart from other fluctuations as implemented in the Monte-Carlo Glauber model. We here emphasize that our present goal is not giving a full theoretical prediction for the observables but quantifying the parton level contributions captured by the CGC model. The flows are characteristic angular distributions defined from the m-particle inclusive spectra, which are in the dipole model given by where B is a dipole model parameter. The dipole description is valid in the large N c limit in which gluons are represented by dipoles. In principle, the dipole distribution should be fixed by the solution of the quantum evolution equation, but in the dipole model the distribution is simply approximated in a Gaussian form. Typically, the Gaussian size B should be of order of the dilute-system (or proton) size ∼ 1 fm ∼ (0.2 GeV) −1 , and we take √ B = 2 GeV −1 . We note that, if (multiple) gluon radiations occur, the above formula would be modified. Such processes (after the collision) are parametrically suppressed by higher powers of JHEP11(2017)114 the strong coupling constant, but may be collinearly enhanced, leading to a parton shower. In the present analysis we implicitly assume that the initial parton correlations are not diminished by the hadronization process. The gluon splitting processes are also possible before the collision, but they should be identified as the parton evolution in the projectile wave-function, which is here approximated by the dipole model.
The general analysis for the flow properties is presented in ref. [48] and the n-th moment of the m-particle correlation is introduced as where φ i represents the azimuthal angle, i.e. p ⊥i = |p ⊥i |(cos φ i , sin φ i ). Then, in the dipole model, we can perform the momentum integrations to find the following expression, D(x ⊥j , y ⊥j ) .

(6.3)
Here, using the regularized generalized hypergeometric function, we defined, where θ is the azimuthal angle of x ⊥ . In particular, we will frequently use the n = 0 function for the normalization, which is given by It is important to notice that ± is irrelevant for n = 0 and there is no angular dependence any more in K (±) 0 (x ⊥ ). Also, we must point out thatδ (2) pmax (x ⊥ ) should approach δ (2) (x ⊥ ) in the p max → ∞ limit.
Let us first consider the case with m = 2 using our large-N c expansion. The N 0 c order term in eq. (4.14) does not contribute to κ n {2} due to the phase factor in K (±) n (x ⊥i − y ⊥i ). As a result, we can write κ n {2} using the N −2 c order results as where .  n . The last part of the integrand is a function of modulus of various combinations of x ⊥1 , x ⊥2 , y ⊥1 , y ⊥2 . Here, it is crucially important to understand that any term in the integrand which is factorized into a function of |x ⊥i − y ⊥i | alone would vanish due to the phase factors in K (±) n (x ⊥i − y ⊥i ) in the factorized integrations. Because there is no finite contribution of disconnected parts in the two particle correlation, we can immediately compute the two particle flow harmonics, v n {2}, from where the denominator is obtained with eq. (6.5), i.e. if we keep using the expanded expression up to the N −2 c order for later convenience, we have Here, we defined in the denominator is in principle irrelevant since it gives a higher order correction which we should neglect.
We summarize our numerical results in figure 3. We have performed the 8 dimensional numerical integration with respect to {x ⊥i , y ⊥i } using the Monte-Carlo method by taking 10 8 sampling points. To draw figure 3 we choseΛ = 0.241 GeV and p max = 2 GeV in accord with ref. [43]. We also note that, only in this section, we change the definition of the saturation momentum from our original Q s in eqs. (2.2) and (2.3) to newQ s defined bȳ

JHEP11(2017)114
according to refs. [42,43]. Since there is no confusion, in this section, we will omit bar and simply denote Q s to meanQ s . Then, we can make a direct comparison of our outputs to figure 1 of ref. [42]. The dashed curved in figure 3 must precisely reproduce figure 1 of ref. [42]. At a glance of our numerical calculations we see quantitatively good agreement. The most interesting question is how useful our large-N c formulae (4.14) can be for the N c = 3 case. The comparison between the dashed curves (full analytical results) and the solid curves (large-N c approximations) in figure 3 concludes that the errors are of only a few (at most ∼ 5) % level except for the n = 3 case. Next, it is intriguing to see what happens for m = 4. In this case, κ n {4} involves K n (x ⊥4 − y ⊥4 ). It is then easy to understand that our formula of the N −2 c order in eq. (4.14) is insufficient to get nonvanishing contributions. Terms of the formula (4.14) are functions of, say, x ⊥1 , x ⊥2 , y ⊥1 , y ⊥2 for i = 2 and j = 1, and then the angle integrations of x ⊥3 − y ⊥3 and x ⊥4 − y ⊥4 become zero. Therefore, one permutation is not enough, but two permutations are necessary to shuffle all x ⊥1 , x ⊥2 , x ⊥3 , and x ⊥4 ; we must go to the next N −4 c order to compute a first nonzero term in κ n {4}.
It is a straightforward but tedious calculation to pick all the N −4 c order terms up from the expansion in eq. (4.11). We can slightly simplify the problem by discarding terms which do not contribute to κ n {4}. The N −4 c order terms generally contain the product of the interaction matrix elements like but we already saw that, if an unexchanged pair exists, the angle integration is vanishing. For example, if the above matrix elements are ∼ V  After long calculations we arrive at a final form which turned out to be factorized as  which starts differing from anticipated eq. (6.18) for m > 4. Thus, to distinguish the connected contribution from the normalization effect, one can test the N c scaling properties as in eqs. (6.18) and (6.21) and also check the p max dependence since δv n {m} from the normalization effect is very sensitive to p max /Q s as perceived from figure 4.

Conclusions
We have established general formulae and machineries to compute 2n-point Wilson line (or n dipole) correlators with the color group representation (N c ⊗N c ) n in the McLerran-Venugopalan model. The color structure accommodates a huge representation but the nonzero contribution to the Wilson line correlators reduces to the n! × n! matrix, whose bases correspond to the color singlets. In particular, we have derived the explicit expression of the matrix elements [see eqs. order for the dipole correlators as given in eq. (4.14). As a check of the validity, we have compared results from the exact answers and those in the large-N c expansion for the two-particle flow harmonics, v n {2} (n = 2, 3, 4, 5), which shows quantitatively good agreement. We have continued our large-N c expansion to higher orders to discuss the flow harmonics with more particles. Then, we have discovered the N c scaling as v n {m} ∼ N −2+2/m c even beyond the glasma graph approximation but in the full MV model. We have also pointed out that a slightly different N c scaling could emerge from the normalization effect at finite cutoff of the transverse momenta of integrated particles.
Although we focused on only the dipole correlators in the present paper, our general formulae also provide us with useful approaches to evaluate Wilson line correlators in channels relevant for the particle production rate in the p-A collision generally. In fact, not only fundamental but also adjoint Wilson line correlators appearing in the multi-gluon production can be derived from our results in eq. (2.1) through the relation [ . As we emphasized, our scheme of the large-N c expansion takes a nice form which is easily implemented in numerical algorithms to go to arbitrarily higher orders. Such higher order numerical evaluations are left as an intriguing future problem. It would be also an important question to think about generalizations beyond the MV model. Further systematic considerations on the Wilson line correlators should deserve more investigations in the future.