Next-to-leading order Balitsky-Kovchegov equation beyond large $N_\mathrm{c}$

We calculate finite-$N_\mathrm{c}$ corrections to the next-to-leading order (NLO) Balitsky-Kovchegov (BK) equation. We find analytical expressions for the necessary correlators of six Wilson lines in terms of the two-point function using the Gaussian approximation. In a suitable basis, the problem reduces from the diagonalization of a six-by-six matrix to the diagonalization of a three-by-three matrix, which can easily be done analytically. We study numerically the effects of these finite-$N_\mathrm{c}$ corrections on the NLO BK equation. In general, we find that the finite-$N_\mathrm{c}$ corrections are smaller than the expected $1/N_\mathrm{c}^2 \sim 10\%$. The corrections may be large for individual correlators, but have less of an influence on the shape of the amplitude as a function of the dipole size. They have an even smaller effect on the evolution speed as a function of rapidity.


I. INTRODUCTION
In hadronic collisions at high energies, large gluon densities are created by the emission of soft gluons carrying a small fraction of the longitudinal momentum of the parent [1]. Nonlinear dynamics of gluons becomes important in such an environment, where parton densities eventually grow to become on the order of the inverse of the QCD coupling α s . To describe QCD in this region, the Color Glass Condensate (CGC) effective field theory [2] has been developed.
In the CGC framework, cross sections for various scattering processes can be expressed in terms of correlators of Wilson lines. A Wilson line describes the eikonal propagation of a parton in the strong color field of the target. The energy dependence of the target color fields, and thus cross sections, is obtained by solving the so-called Jalilian-Marian-Iancu-McLerran-Weigert-Leonidov-Kovner (JIMWLK) equation [3][4][5][6]. This is a perturbative evolution equation that describes the Bjorken-x dependence of a Wilson line. In phenomenological applications, it is usually convenient to work directly in terms of the Wilson line correlators, and to solve instead the Balitsky-Kovchegov (BK) equation [7,8] for the dipole operator (correlator of two Wilson lines), which can be obtained from the JIMWLK equation in the large-N c limit.
The CGC framework has been used extensively in phenomenological applications at leading order (LO) in α s , with the evolution equations resumming contributions ∼ α s ln 1/x to all orders. Running coupling effects derived in Refs. [9][10][11][12] (see also [13]) can also be taken into account. The non-perturbative initial condition for the small-x evolution is obtained by performing fits to the HERA structure function data [1,14], for example in Refs. [15][16][17][18] (see also [19,20]). The obtained initial condition can then be used for various calculations, for example particle production in proton-nucleus collisions [17,[21][22][23][24][25][26][27]. In the future, the nuclear deep inelastic scattering (DIS) experiments at the Electron Ion Collider (EIC) [28,29] in the US, at the LHeC [30] at CERN and at the EicC in China [31] will provide a vast amount of precise data from clean DIS processes. These experiments will be able to probe the nuclear structure where nonlinearities are enhanced by roughly A 1/3 higher densities compared to the proton. Before the EIC, similar studies limited to the photoproduction region can be performed in ultra-peripheral heavy-ion collisions [32,33].
The BK equation is usually solved in the large-N c limit. In the LO case, the large-N c limit makes it possible to express the four-point correlator of fundamental representation Wilson lines in terms of the two-point function.
In detailed numerical studies, it has been shown that the finite-N c corrections are smaller than the naive expectation of O(1/N 2 c ) [52,53]. At NLO, the equation involves six-point functions of fundamental Wilson lines that must similarly be expressed in terms of the two-point function in order to close the equation. The purpose of this work is to see if the finite-N c corrections are similarly small in the case of the NLO equation, where all corrections of the order α 2 s are taken into account. In order to numerically solve the BK equation at finite N c , we use the Gaussian approximation [54][55][56] to derive analytical parametric equations for the six-point correlators in terms of the two-point correlators. We study numerically the finite-N c corrections to the these correlators, and also their effect on the NLO BK evolution. In addition to the BK equation, higher-point correlators are needed in the calculations of multi-particle correlations in the CGC framework, see eg. Refs. [57][58][59][60].
The structure of the paper is as follows. In Section II, we introduce the NLO BK equation and provide both the arXiv:2007.00751v1 [hep-ph] 1 Jul 2020 large-N c and finite-N c expressions for the correlators that will be studied. In Section III, we introduce the Gaussian approximation, explain the diagrammatic notation used in the rest of the paper and then explain the analytical calculation done for finding the parametric equations for the six-point correlators. Section IV contains the numerical results obtained from using the analytical expressions for the six-point correlators to solve the BK equation. Finally, we end with a few concluding remarks and a summary of our main results.

II. THE BK EQUATION AT NLO
For any product of n/2 pairs of fundamental Wilson lines UU † , we use the notation The NLO BK equation in the case of zero active quark flavors (n f = 0) reads [34] where the brackets refer to the expectation value over target color field configurations. The kernels are The convolutions ⊗ in Eq. (2) denote integrations over the transverse coordinate z (in K BC 1 ) or z and z (in K 2,1 and K 2,2 ). We use the notation We note that the kernel proportional to n f is also available [34]. Since the purpose of this work is to study the importance of the finite-N c corrections in the NLO BK equation, we do not include contributions proportional to n f . The finite-N c effects could be expected to be similar to the n f = 0 case.
The Wilson line operators appearing in Eq. (2) are Although the original NLO BK equation in the form presented in Ref. [34] does not contain the subtraction z → z in D 2,2 , we have introduced the subtraction to improve numerical stability. This subtraction term has no effect on the final evolution because the integral of K 2,2 over z vanishes if the Wilson line operator term does not depend on z (see Ref. [34]). We will refer throughout this work to two pieces of the right side of Eq. (2) as the In other words, we separate the terms in the NLO BK equation by the types of Wilson line correlators, not by the order in α s . Thus, the LO-like contribution also includes a significant α 2 s correction. The interpretation of the NLO BK equation is that one considers all possible ways to emit either one or two gluons, at transverse coordinates z and z , from the boosted dipole consisting of quarks at transverse coordinates x and y. The effect of the boost is that instead of the original dipole projectile, the original quarks and the emitted gluons scatter off the target color field. As such, the evolution can be seen to describe the evolution of the projectile probing the target structure. On the other hand, the emitted gluons can also be taken to be a part of the target wave function, in which case the boost corresponds to the evolving target color field as probed by the original projectile. For a more detailed discussion on the NLO evolution in the projectile or target wave function, the reader is referred to Ref. [61].
The NLO BK equation is known to be unstable [62] due to the large contributions enhanced by the large double transverse logarithm ln X 2 /r 2 ln Y 2 /r 2 . We resum these contributions to all orders following the procedure developed in Ref. [63], which was numerically confirmed in Ref. [64] to result in a stable evolution (see also Ref. [65] for an equivalent resummation of the same double logarithms). In addition, we include the running of the QCD coupling by noticing that the terms proportional to the beta function coefficient β in Eq. (3) should be resummed into the running coupling. We implement this resummation by following the Balitsky prescription from Ref. [12]. Both running coupling and double transverse logarithm resummations are included by modifying the kernel K BC The double log corrections to all orders are taken into account by the factor whereᾱ s = α s N c /π. The double logarithm here is x = ln X 2 /r 2 ln Y 2 /r 2 . If ln X 2 /r 2 ln Y 2 /r 2 < 0, then an absolute value is used and the Bessel function is changed from J 1 → I 1 (see Ref. [63]). The scale of the coupling in K DLA is determined by the smallest dipole min{r 2 , X 2 , Y 2 }.
In addition to the double log contributions, one can also resum a set of higher-order contributions enhanced by single transverse logarithms, as shown in Ref. [66]. For the purposes of this paper, this resummation is not necessary and is excluded for simplicity. In this running coupling prescription, we keep the order α 2 s terms in the kernel K BC 1 that are not proportional to the beta function. These are included in the term K fin 1 , which reads The strong coupling constant in the transverse coordinate space is evaluated as where β = (11N c − 2n f )/N c . We use * C 2 = 1 and µ 0 /Λ QCD = 2.5 in our numerical calculations, which freezes the coupling at α s (r → ∞) = 0.762 in the infrared, and c = 0.2 which controls the transition to the infrared region. The initial condition for the BK equation is taken from the McLerran-Venugopalan (MV) model [67,68]. In the MV model, the color charge density is assumed to be a random Gaussian variable, with a zero expectation value and a variance proportional to the local saturation scale Q 2 s . The dipole correlator in the MV model is written as Here, the constant e acts as an infrared regulator. We use Λ QCD = 0.241 GeV and Q 2 s0 = 1 GeV 2 in the numerical analysis. In analytical studies of the correlators of Wilson lines in specific "line" coordinate configurations in Sec. IV A, we use the GBW [69] form for the dipole correlator In principle, the resummation procedure for the double transverse logs would also change the initial condition, as discussed in Refs. [63,64]. However, since the initial condition is a non-perturbative input for the evolution, we consider Eq. (13) to be the non-perturbative initial condition for the resummed evolution as well. For the purposes of this paper, the actual form of the initial condition is not relevant. We compare the finite-N c version of the BK equation presented above to the equation obtained in the large-N c limit, which has been studied numerically in Refs. [62,64]. In this limit, one can drop operators suppressed by 1/N c and the correlators in Eq. (2) become with [9,13] would be C 2 = e −2γ E ≈ 0.32. We use a larger value C 2 = 1 which results in slightly slower evolution, as the C 2 is usually taken to be a free parameter controlling the coordinate space scale.

A. The Gaussian approximation
At large N c , all Wilson line operators present in the NLO BK equation can be expressed solely in terms of dipole correlators, as can be seen from Eqs. (15) and (16). At finite N c , on the other hand, the higher-point functions S (2) x,z S and S (6) x,z,z ,y,z,z are needed. An expression for the four-point function in terms of two-point functions has been derived using the Gaussian approximation (see e.g. [55]). This makes it possible to obtain a closed form for the LO BK equation at finite N c . The purpose of this work is to compute also the six-point functions using the Gaussian approximation, in order to obtain a closed finite-N c BK equation at NLO accuracy.
In the Gaussian approximation, all correlators are parametrized by a single two-point function, and all higher-point functions can then be expressed in terms of this function. The initial condition for the small-x evolution is usually assumed to be Gaussian (e.g. as in the MV model), but it is not clear a priori that the Gaussian approximation is valid after the evolution. However, numerical studies of the JIMWLK equation [70,71] have not found any indication of major effects breaking the validity of this approximation. We use the diagrammatic notation of Refs. [53,55,72], in which Wilson lines are denoted as The projectile transverse coordinate is x, the lightcone time axis runs from right to left and the blue vertical line represents the target background field. In the Gaussian approximation, the correlator for some Wilson line operator O[U ] is approximated as an integral over a parametrization rapidity η of a single two-point correlator G u1,u2 : The transverse integrals are denoted as u = d 2 u and L a u is a Lie derivative that acts on Wilson lines according The structure as an exponential of a two-point function (as in the MV model [67]) is what makes this a "Gaussian" approximation. In practice, for gauge invariant (color singlet) operators, the two-point function G u1,u2 always appears in the linear combination For the integrand, we use the notation G := ∂ η G. Physical observables only depend on the integrated G and not on the integrand G . Thus, for the purpose of our calculation, where we need to relate higher-point functions of Wilson lines to the two-point function, there is some freedom in choosing the parametrization rapidity. We use this freedom in such a way that the η and transverse coordinate dependences of G factorize * , as is usually done when employing the Gaussian approximation (see e.g. [54,56,58,70,[74][75][76]) in phenomenological applications. We henceforth omit the explicit η dependence of G for brevity.
As an example of the procedure for finding the parametric equation for a correlator using Eq. (21), we illustrate the steps for the dipole operator S (2) x,y . A more practical form of Eq. (21) is to write the integral over the parametrization rapidity η in differential form. Then, we have where M is the so-called transition matrix formed by operating with the argument of the exponential in Eq. (21) on the operator of interest. In this case, there are four contributions from acting with the Lie derivatives on a * This assumption has the effect that the transition matrices M(η) (introduced below in Eq. (30)) at different rapidities η commute with each other. This turns the path ordered exponential of M(η) into a normal exponential. This, in turn, makes it possible to relate higher-point functions to the two-point function without any further assumptions about the η dependence of G . In the terminology of Ref. [55], we use a "rigid exponentiation" instead of the "Gaussian Truncation". The "Gaussian Truncation" would imply equating the parametrization rapidity η with the evolution rapidity Y . See the related discussion in Refs. [72,73].
product of two Wilson lines: From this sum, we need to factorize out the original operator U x ⊗ U † y , so we use the Fierz identity There is only one way to join the endpoints of the Wilson lines into a singlet operator. This is to trace over them, i.e. wedging them between 1 2Nc . Since the operators and were normalized, the initial condition at η = 0 for the differential equation (25) is given by trivial Wilson lines equal to the identity matrix in the absence of a color field. This is the well-known parametric equation for the dipole correlator [55] in the Gaussian approximation.
In the case of n-point correlators larger than the dipole, the operator O[U ] in Eq. (21) is actually an n×n matrix of correlators, denoted A(η), and Eq. (25) becomes an n × n matrix differential equation By construction, M is a symmetric matrix, so there are at most n i=1 i = n(n + 1)/2 distinct elements, not n × n. For example, a product of six Wilson lines is repre-sented in this notation as (the haphazard assignment of coordinate labels is convenient for the NLO BK equation and will become clear when constructing the transition matrix). The notation here means that this product is actually a matrix with six open indices on the left and another six on the right; we denote them as Since there are six possible ways to join the endpoints of these Wilson lines to form a correlator, Eq. (30) is a six-by-six matrix differential equation, as opposed to the much simpler one-dimensional problem illustrated for the dipole correlator.
In analogy to the procedure for the dipole operator, the procedure to use Eq. (30) to find parametric equations for the six-point correlators is as follows: and the corresponding element for the other end of the Wilson line is then The prefactor is a normalization constant found by squaring the basis element.

Construct the correlator matrix
3. Construct the transition matrix M by summing (for each element in A) all possible one-gluon diagrams obtained with the double Lie derivative operator and rewriting the result in terms of elements of A.
4. Solve Eq. (30) by exponentiating M to find expressions for each element in A, using as an initial condition the correlator matrix A corresponding to Wilson lines equal to the identity matrix.

B. Choosing a basis
Starting from a product of six Wilson lines, there are six ways to form multiplets by joining endpoints in all possible ways: .
The simplest way to construct an orthonormal basis from these would be to use color algebra to choose The Nc . The next step would be to use this basis to construct the correlator matrix A and the transition matrix M. However, doing so results in a matrix differential equation ∂ η A(η) = −M(η)A(η), whose complicated solution is the matrix exponential of a six-by-six matrix M.
For our case, a better way to proceed is to exploit the structure of the six-point correlators that are actually needed for the NLO BK equation. Since there are only four distinct coordinates in these particular correlators, we make the coordinate assignments It is easy to see from this that there is one way to join the endpoints such that in the limit v → z , w → z, four Wilson lines cancel (due to unitarity). The result simplifies to a single trace: So choosing 1 √ N 3 c as one of our basis elements allows one dimension of our six-dimensional space of operators to decouple, giving the equation for the dipole correlator. Similarly, the choice of two more particular basis elements results in two correlators that reduce to fourpoint functions; one due to the limit v → z and the other due to the limit w → z. Thus, we can expect to choose a further two basis elements such that two more dimensions decouple from the remaining five, corresponding to the equation for the four-point operators. These two basis elements can be chosen as We will choose the remaining three basis elements such that they are orthonormal to the three already chosen, resulting in the final basis vector Since this basis is orthonormal, the correlator matrix at the initial condition A(η = 0) will just be the identity matrix.

C. Constructing the correlator matrix and the transition matrix
Due to this choice of basisB, the full matrix differential equation (30) now decouples into three independent equations. This allows us to forego exponentiating a sixby-six matrix; at most we will need to exponentiate a three-by-three matrix, which can be done analytically.
To form the correlator matrix A, we take the product To form the transition matrix, we act with the argument of the exponential in Eq. (21) on the then wedge the result between the basis vectors and set w → z and v → z : Diagrammatically, this is equivalent to summing all possible ways of attaching one gluon line on the operator , using the Fierz identity to replace the gluon vertices, and finally closing the Wilson line endpoints on the left and right using the basis vector. For example, the element (6, 6) of the correlator matrix A is We then sum all the diagrams in which a gluon is attached to this diagram so that it joins any two of the six Wilson lines on the right of the target interaction. Using the Fierz identity in Eq. (28), we may write the resulting expression in terms of diagrams with no gluons. After making the substitutions w → z and v → z , the result will be a linear combination of elements of the operator matrix A(η).
where the subscripts refer to the dimension of the submatrix.
The first (one-dimensional) transition sub-matrix is and upon exponentiation, gives the parametric equation for the dipole correlator as shown in Eq. (29). When inverted, this equation can be used to express the twopoint function G x,y in terms of the dipole correlator. This will be needed to evaluate the higher-point functions in terms of the dipole S x,y . The second (two-dimensional) transition sub-matrix is where The matrix differential equation then gives a coupled system of 2×2 differential equations, out of which 2 are linearly independent, corresponding to the fact that the same transition matrix operates separately on each of the two columns of A 2 . The exponential solution for this system of equations gives the known parametrization for the four-point correlator with one repeated coordinate [55] S (2) x,z S (2) z,y = 1 N 2 c e −CFGx,y + 2C F N c e −CFGx,y e − Nc 2 (Gx,z+Gy,z−Gx,y) . (48) The third and final (three-dimensional) transition submatrix is where and the primes on the Γ's in Eq. (49) denote derivatives in η. Exponentiating this matrix is the last step required to get expressions for the remaining six-point correlators in A 3 .

D. Exponentiating the transition matrix M3
In order to obtain the six-point functions, it is necessary to solve the differential equation Solving this equation is equivalent to exponentiating the matrix M 3 , as shown above in the cases of the two-and four-point functions. To exponentiate M 3 , we consider two different cases: Γ 2 = 0 and Γ 2 = 0. When Γ 2 = 0, M 3 in Eq. (49) becomes diagonal and we directly obtain When Γ 2 = 0, the matrix elements of A 3 are calculated by matrix-exponentiating the full M 3 in Eq. (49), giving Here, z i are the roots of the cubic polynomial They are where Notice that c 3 may be complex. The functions of the roots that appear in Eq. (55) are Despite the fact that some of these expressions are complex, the final expressions for each element in A 3 are in fact real, as they should be.

E. Extracting six-point correlators needed for NLO BK
Eq. (55) gives the analytical expressions for the correlators formed using basisB, solely in terms of the parameter G (which can be used to relate these expressions to the dipole via Eq. (29)). For example, However, the two correlators required in Eqs. (7) and (8) are These are not explicitly any of the elements of matrix A 3 , since and are not basis elements inB. Instead, they are linear combinations of the elementsB i contained inB: Using these two expressions, it is simple to get the final equations for the two six-point correlators needed. In terms of the elements of the correlator matrix A 3 given in Eq. (55) (see Appendix A for detailed expressions) they are x,z,z ,y,z,z and S (6) x,z,z ,y,z,z = − S (2) x,y + S (2) x,z S (2) z,y + S (2) x,z S . (76) Equations (75) and (76) are the final two expressions needed to solve the NLO BK equation at finite N c ; they are the main analytical results of this work. It is now possible to express these six-point functions entirely in terms of dipole correlators using Eq. (29). This makes it possible to write the NLO BK equation from Eq. (2) solely in terms of dipole correlators. In such a closed form, it can be solved directly, as was done in the large-N c case in Refs. [62,64].
To verify the validity of Eqs. (75) and (76), we perform two checks. Firstly, the Gaussian approximation has the built-in property that it should be consistent in color algebra. This means that taking any coincidence limit in which coordinates are made equal in Eqs. (75) and (76), should reduce them to the relevant expressions for the lower-point functions Eqs. (29) and (48). For example, setting z → x and z → y in Eq. (75), we reproduce the equation for the dipole (29), as expected.
Secondly, when Eq. (75) is taken in the dilute limit, where the Wilson lines are expanded as Eq. (75) should be the same up to order λ 2 as the parametric equation for the large-N c counterpart operator.
In the case of correlator S (2) x,z S z,z S (2) z ,y , the large-N c result is just the factorized product of dipole correlators z ,y = e −CF(Gx,z+G z,z +G z ,y ) . (79) After some algebra, Eq. (75) can be shown to give the same result up to order G.
We note that in Ref. [59], correlators of up to eight Wilson lines have been calculated at finite N c . The difference between that work and ours is that the authors there are solving the system for a general configuration of coordinates, where it is difficult to find a basis such that the transition matrix would become block diagonal. Consequently, an analytical approach as presented in this paper is not possible. Instead, the authors numerically exponentiate the transition matrix, which is a much more expensive computational procedure than what is needed here.

IV. NUMERICAL RESULTS
We now study numerically the obtained six-point correlators, Eqs. (75) and (76). In particular, we are inter-ested in the effects of the 1/N 2 c suppressed contributions included in these six-point correlators, compared to the large-N c version in Eq. (16), which was used previously in numerical studies of the NLO BK equation. We will first study these operators in a specific coordinate configuration (with the GBW parametrization for the dipole). We will then integrate the operators over the gluon coordinates z, z and study the BK evolution starting from an MV model initial condition.

A. Correlators in a line configuration of coordinates
As a baseline for comparison of the six-point correlators in the NLO BK integrand, we consider first the four-point correlator S (2) x,z S (2) z,y that appears in the LO BK equation. We compare the full finite-N c result to its large-N c limit S (2) x,z S (2) z,y . The finite-N c correlator is evaluated by applying Eq. (48). For the dipole operator S (2) , we use the GBW form given in Eq. (14). We consider the following, specific but not atypical, configuration of coordinates: x, y and z in a line along the horizontal axis of transverse coordinate space, as shown in Fig. 1. The distance between points y and z is denoted by a, the distance between points x and z is 2a, and the distance between x and y is 3a. We have confirmed that the rough relative magnitude of the finite-N c effects of the results shown in this subsection are not specific to the actual chosen geometric configuration. To show the results as a function of the dimensionless distance scale aQ s , we define the saturation scale Q s as S (2) x,y (x−y) 2 =2/Q 2 In Fig. 2, the four-point correlator is shown both at finite N c and at large N c as a function of the distance aQ s . Additionally, the magnitude of the finite-N c correction is shown as a difference between the finite-N c and large-N c results, denoted by S (2) S (2) − S (2) S (2) . The finite-N c correction to the four-point correlator S (2) S (2) is found to be negligible. At typical aQ s = 1, the relative finite-N c correction ( S (2) S (2) − S (2) S (2) )/( S (2) S (2) ) is approximately 5%. The relative correction becomes more important at large aQ s , in the region which gives only a negligible contribution to the BK evolution. We will return to the discussion of the finite-N c corrections at aQ s 1 later, when evaluating the six-point functions. Also shown in Fig. 2 is the full LO-like operator factor D 1 (see Eq. (6)) from the BK equation, both at large and finite N c . The difference between the finite-N c and large-N c results is the same as the difference for the four-point correlators. The fact that the finite-N c corrections are smaller than ∼ 1/N 2 c is not surprising, as these corrections to the LO BK equation are known to be small [53].
Next, we choose for the four coordinates present in the NLO-like operators in the BK equation a similar line con-  Fig. 2, it is still negligible at aQ s 1. On the other hand, the finite-N c corrections clearly dominate in the region aQ s 1, the relative contribution from 1/N 2 c suppressed terms being approximately 40% at aQ s = 1. A similar, although numerically smaller, effect was observed in the four-point function studied above. This can be understood as follows. When 2a 1/Q s , the color fields at points x and z, as well as at y and z are uncorrelated. However, at finite N c , the six-point function is also sensitive to the color field correlations between points x and z , as well as between y and z that belong to different dipoles and are thus not correlated in the large-N c limit. When |x − z | = |y − z| 1/Q s , these correlations do not vanish and actually dominate the full six-point function.
Also shown in Fig. 4 for comparison is the other 1/N 2 c suppressed six-point correlator 1 N 2 c S (6) present in the NLO BK integrand at finite-N c (cf. Eq. (7)). This is plotted using Eq. (76). We can see that the contribution of the six-point function S (6) /N 2 c to D 2,2 is similar in magnitude as that of the finite-N c corrections to the dipole cubed operator S (2) S (2) S (2) .
In Fig. 5, we use all the above mentioned correlators to plot the NLO-like factors D 2,1 and D 2,2 , as defined by Eqs. (7) and (8), respectively. Since both quantities reduce to the same expression in the large-N c limit, only one curve is shown for the large-N c case. The dashed curves show the differences, i.e. the finite-N c corrections to D 2,1 and D 2,2 . As already seen when studying the six-point correlators, the finite-N c corrections are negligible at aQ s 1, but become numerically important when aQ s 1. In comparison to the dashed curves in Fig. 2 for the LO-like case, we see that the finite-N c corrections in the NLO-like case are larger. At aQ s = 1, the finite-N c corrections to the operators D 2,1 and D 2,2 are approximately 20% and 16%, respectively. In comparison, the LO-like operator D 1 shown in Fig. 2 has a finite-N c correction of approximately 8%. When considering the full NLO BK evolution, one should keep in mind that the evolution is driven by the dipole sizes r 1/Q s . As such, even though the finite-N c corrections can be large at r = 1/Q s , the actual effect of the 1/N 2 c suppressed contributions to the small-x evolution can be smaller. The NLO BK evolution at finite-N c is studied in the next section.

B. BK evolution at finite Nc
The line configuration studied in the preceding subsection is a typical one, and the curves shown represent the expected behaviour of the correlators. Equipped with this, we move on to studying the full BK equation using the MV model initial condition, shown in Eq. (13), with Q 2 s,0 = 1 GeV 2 . Since the finite-N c corrections to the individual operators have been found to be small (except at large distances which do not significantly contribute to the BK evolution), we expect the finite-N c corrections to remain small when performing integrations over gluon coordinates z and z in Eq. (2). In Fig. 6, we show the relative evolution speed 1 N ∂ Y N , where the dipole amplitude N = N x,y is defined as x,y . This is obtained by integrating the full right side of Eq. (2), first using the large-N c expressions for the correlators, then again using the finite-N c expressions. These are shown separately for the LO-like contribution (only the term containing D 1 in the integrand in Eq. (2)) and the NLO-like contribution (only the terms containing D 2,1 and D 2,2 ). As expected from the line configuration studies, we see that the finite-N c corrections for the NLO-like terms are slightly larger, but of the same order of magnitude as the finite-N c corrections for the LO-like terms. The finite-N c corrections vanish when the parent dipole size r is small, and are most important at rQ s ∼ 1, as expected from the line configuration analysis presented above.
In Fig. 7, we plot the difference between the large-N c and finite-N c cases, separately for the LO-like and NLO-like terms. This shows more clearly that the difference for the NLO-like terms is of the same order of magnitude as for the LO-like terms. We also note that the difference has the opposite sign in the LO-like and the NLO-like terms. Consequently, a part of the difference cancels in the total evolution speed. At rQ s = 1, the relative finite-N c correction is approximately 8% in the LO-like contribution and 13% in the NLO-like con- tribution. The relative magnitude of the total 1/N 2 c suppressed contribution is 5%, which is somewhat smaller than the expected correction of 1/N 2 c ∼ 10%. Finally, the last thing left to study is to move beyond the initial condition and determine how the finite-N c corrections behave under the NLO BK evolution. In Fig. 8, we show the ratio of the dipole amplitudes N obtained by solving the full NLO BK equation at finite N c to that at large N c . At r 1/Q s , when the details of the initial condition are lost and one enters the geometric scaling region, the difference between the large-N c and finite-N c cases evolves only very slowly. At small dipoles, the ratio grows approximately linearly in Y . The fact that the total finite-N c correction is positive at small dipole sizes and negative at large dipoles, as seen in Fig. 7, is found to hold also asymptotically after many units of rapidity evolution.
The evolution speed of the saturation scale, ∂ Y ln Q 2 s , is shown in Fig. 9. Similarly to what is seen in the dipole amplitude plot in Figure 8, we see from this figure that the finite-N c corrections are more important at the initial condition, slowing down the evolution of Q 2 s by approximately 5%. Later in the evolution, where the solution approaches the asymptotic shape of the BK evolved dipole, the difference becomes smaller -of the order of 1%. Consequently, even at the initial condition (and especially when the details of the initial condition are lost) the finite-N c corrections to the evolution speed of Q s are found to be significantly smaller than the naive expectation of 1/N 2 c at NLO.

V. CONCLUSIONS
In this work, we have studied the six-point correlators in the NLO BK equation using the Gaussian approximation. This allowed us to express these higher-point correlators in terms of the dipole operator. In using our analytical results, we have seen numerically that the overall finite-N c corrections to the NLO-like part of the BK equation are somewhat smaller than what is naively expected. However, one needs to state the actual quantity being compared in order to quantify this correction.
When correlators are considered between Wilson lines separated by large distances relative to 1/Q s , 1/N 2 c suppressed corrections may be considerable. Despite these potentially large corrections to individual correlators, these configurations do not contribute much to the right side of the BK equation. Therefore, we find a somewhat smaller, although still significant, effect on the shape of the dipole amplitude as a function of r. The finite-N c corrections are watered down further when one considers the evolution speed of Q s as a function of rapidity, especially once the evolution settles towards its asymptotic form away from the initial condition. In general, finite-N c corrections need to be considered carefully when evaluating the NLO BK equation, since they may have a non-negligible effect at the required accuracy. The content of this article does not reflect the official opinion of the European Union and responsibility for the information and views expressed therein lies entirely with the authors. Computing resources from CSC -IT Center for Science in Espoo, Finland and from the Finnish Grid and Cloud Infrastructure (persistent identifier urn:nbn:fi:research-infras-2016072533) were used in this work.

A. CORRELATOR MATRIX A(η)
We give here the explicit expressions for the operators contained in the correlator matrix A(η). Since the transition matrix M block-diagonalizes in basisB, as explained in Section III C, we are only interested in the corresponding block-diagonalized matrix The one-dimensional sub-matrix is The two-dimensional sub-matrix is where A (1,1) 2 The three-dimensional sub-matrix is