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/N_c$ where $N_c$ is the color number. As a phenomenological application we calculate the flow harmonics $v_n\{m\}$ in the dipole model and discuss the $N_c$ scaling.


Introduction
Quantum Chromodynamics (QCD) has a special property called the asymptotic freedom that implies that the strong coupling constant runs to a smaller value in parton reactions involving harder momentum scales. In this way the perturbative calculation of QCD is a reliable theoretical description in high-energy nuclear physics. In reality, however, a naïve perturbative expansion breaks down for diffractive type reactions in which exchanged momenta are not necessary hard as compared to the collision energy. Then, it is indispensable to make a resummation over large logarithmic enhancement factors which appear kinematically. Such a resummation leads to a picture of increasing parton or gluon density with increasing scattering energy, and, one would anticipate that the parton density would eventually enter a new regime where the effect of parton overlapping is significant. In such a regime of dense partons, QCD is still perturbative in a sense that the running coupling constant g is small, but is highly non-linear due to the gluon amplitude A µ as large as ∼ 1/g [and thus gA µ ∼ O (1)]. This non-linear dynamics is a manifestation of the parton or gluon saturation (see Refs. [1][2][3][4][5] for pioneering extensions including the non-linearity, see also Refs. [6,7] for recent reviews).
The virtue of the gluon saturation is that physical observables exhibit universal behavior in terms of scaling variables, from which the saturation momentum, Q s , can be experimentally fixed [8]. The theoretical framework in the saturation regime to compute physical observables as functions of Q s has been well established and known as the Color Glass Condensate (CGC) [9][10][11]. The CGC effective theory is elegantly formulated in the renormalization group language [12][13][14], in which the soft gluons are given by the classical Yang-Mills fields from the hard parton color source [15] whose transverse density is 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, 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 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 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 Sec. 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 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 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, 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; , so that we can again give a schematic representation as Also, another useful formula is 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 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 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.
The definition of the dipole operator is 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 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 Fig. 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 Fig. 1.
We next see the Q s dependence of the n = 2 dipole correlator together with the Abelian approximation, as depicted in Fig. 2. We introduce the Abelian approximation as employed in Ref. [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 extrapolation from the n = 1 answer. Figure 2 clearly shows that a small disagreement magnified in Fig. 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}, according to Refs. [42,43]. 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, which is typically of the order of the nucleon size ∼ 1 fm ∼ (0.2 GeV) −1 , and we take √ B = 2 GeV −1 . 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, 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 where . (6.7) One could think of D 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  . The saturation momentum Q s is the one defined in Eq. (6.11) as in Ref. [43].
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 pmax (x ⊥ − y ⊥ ) . in the denominator is in principle irrelevant since it gives a higher order correction which we should neglect.
We summarize our numerical results in Fig. 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 Fig. 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ȳ 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 Fig. 1 of Ref. [42]. The dashed curved in Fig. 3 must precisely reproduce Fig. 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 Fig. 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  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.